Thomas Cucchi
email: thomas.cucchi@mnhn.fr

Barbara Stopp

Renate Schafberg

Joséphine Lesur

Alexandre Hassanin

Jörg Schibler

Taxonomic and phylogenetic signals in bovini cheek teeth: Towards new biosystematic markers to explore the history of wild and domestic cattle

Keywords: Bovini, domestication, zooarchaeology, geometric morphometrics, teeth, phylogenetic signal, taxonomy

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

The domestication of aurochsen in several parts of Eurasia had a dramatic effect on human history. Their domestic forms became the symbols of prestige, leadership and wealth, while their strength intensified agricultural production, driving the rise of cities and civilisations [START_REF] Fagan | The intimate bond: how animals shaped human history[END_REF]. Nowadays, they provide nearly 30% of the world's meat and over 87% of the world's milk [START_REF] Scherf | World watch list for domestic animal diversity 3rd edition ed., Food and Agriculture Organization of the United Nations (FAO)[END_REF]. They are also one of the main sources of unsustainability in our food production system, since 26% of the ice-free land on our planet is used for grazing livestock and a third of cropland is used for livestock feed production (www.fao.org/nr/sustainability/sustainability-and-livestock). Aurochs domestication also played a major role in the evolution of the tribe Bovini (Figure 1). Two domestic forms stemmed from aurochs: common cattle (Bos primigenius taurus) and zebu cattle (Bos primigenius indicus), one from the yak: the domestic yak (Bos mutus grunniens), one from the gaur: the gayal (Bos gaurus frontalis), and one from the Asian water buffalo: the river buffalo (Bubalus arnee bubalis) [START_REF] Hassanin | Systematic and evolution of Bovini[END_REF]. Currently, due to habitat loss, these wild ancestors, if not extinct or endangered, have been reduced to small pockets of populations in confined areas (Figure 1). The worldwide dispersal of common and zebu cattle induced adaptations to highly diverse environments, such as lowlands and highlands or humid and arid conditions, driving the emergence of different phenotypes (e.g. [START_REF] Felius | On the Breeds of Cattle -Historic and Current Classifications[END_REF]. This variety of characteristics evolved over thousands of years, but the emergence and development of welldefined, genetically isolated breeds was mainly established from the Industrial Revolution onwards [START_REF] Ajmone-Marsan | On the Origin of Cattle: How Aurochs Became Cattle and Colonized the World[END_REF].

So far, the domestication of aurochs is the most documented history from archaeological and genetic evidence. The earliest evidence (10,500 BP) of the human control of aurochsen based on male body size reduction comes from the northern Fertile Crescent [START_REF] Helmer | Identifying early domestic cattle from Pre-Pottery Neolithic sites on the Euphrates using sexual dimorphism[END_REF], Peters et al 2005). Since common or taurine cattle have a potential progenitor in the form of European aurochs, it is still debated how much local European aurochs have contributed to the genetic make-up of modern cattle (e.g. [START_REF] Cymbron | Microsatellite diversity suggests different histories for Mediterranean and Northern European cattle populations[END_REF][START_REF] Beja-Pereira | The origin of European cattle: Evidence from modern and ancient DNA[END_REF][START_REF] Tresset | Early diffusion of domestic bovids in Europe: an indicator for human contacts, exchanges and migrations?[END_REF][START_REF] Schibler | Incorporation of aurochs into a cattle herd in Neolithic Europe: single event or breeding?[END_REF][START_REF] Scheu | The genetic prehistory of domesticated cattle from their origin to the spread across Europe[END_REF][START_REF] Upadhyay | Genetic origin, admixture and population history of aurochs (Bos primigenius) and primitive European cattle[END_REF]. The second domestication centre led to the origin of zebu cattle, also referred to as humped cattle. A general consensus stemming from archaeological [START_REF] Meadow | The Origins and Spread of Agriculture and Pastoralism in Northwestern South Asia[END_REF] and genetic evidence (e.g. [START_REF] Bradley | Genetics and the Origins of Domestic Cattle[END_REF] considers the Indus Valley as the most likely centre of domestication for this animal. However, other regions in the Indian subcontinent are also considered potential domestication centres [START_REF] Chen | Zebu Cattle Are an Exclusive Legacy of the South Asia Neolithic[END_REF]. The third domestication centre for aurochs has been claimed as Africa [START_REF] Pérez-Pardal | Multiple paternal origins of domestic cattle revealed by Y-specific 29 interspersed multilocus microsatellites[END_REF], but genetic evidence suggests that the first domestic cattle here were actually taurine cattle from the Middle East, with an admixture from local African aurochs [START_REF] Pitt | Domestication of cattle: Two or three events?[END_REF]. In mainland and island south-east Asia, the domestication history of the yak in the Tibetan Plateau or the Asian water buffalo, the gaur and the banteng, mainly relies on genetic and genomic investigations from current wild and domestic populations [START_REF] Bradley | Genetics and the Origins of Domestic Cattle[END_REF]). Although genetic and paleogenetic studies have provided huge insights into the broad origin and trajectories of taurine and zebu cattle, a refined understanding of the timing and mechanisms of their emergence remains to be undertaken in south-western Asia [START_REF] Arbuckle | Documenting the initial appearance of domestic cattle in the Eastern Fertile Crescent (northern Iraq and western Iran)[END_REF]. In southeast Asia, studies on the domestication process of wild cattle through archaeological evidence is progressing, with the first ancient mtDNA sequences just being published [START_REF] Siripan | Origin of prehistoric cattle excavated from four archaeological sites in central and 30 northeastern Thailand[END_REF].

Along with genetic markers, documenting the domestication pathway of wild cattle could greatly benefit from more powerful morphometric tools to document, from the fossil record, the distribution range of wild cattle species and sub-species before the domestication process was initiated, and identify intraspecific morphological changes induced by the domestication process. The morphological markers used to identify cattle domestication from aurochs and subsequent livestock diffusion have so far relied on body-size reduction (Peters et al. 2005;[START_REF] Helmer | Identifying early domestic cattle from Pre-Pottery Neolithic sites on the Euphrates using sexual dimorphism[END_REF]2015[START_REF] Arbuckle | Data Sharing Reveals Complexity in the Westward Spread of Domestic Animals across Neolithic Turkey[END_REF]. This body-size reduction, used as domestication syndrome [START_REF] Hemmer | Domestication. The decline of Environmental perception[END_REF], is a pleiotropic consequence of the selection for behaviours facilitating adaption to the human environment [START_REF] Trut | Animal evolution during domestication: the domesticated fox as a model[END_REF][START_REF] Wilkins | The "Domestication Syndrome" in Mammals: A Unified Explanation Based on Neural Crest Cell Behavior and Genetics[END_REF]. However, difficulty in disentangling human intentionality from other potential biotic and abiotic factors [START_REF] Maurer | THE MICRO AND MACRO IN BODY SIZE EVOLUTION[END_REF] makes it a poor phenotypic marker of early domestication [START_REF] Helmer | Identifying early domestic cattle from Pre-Pottery Neolithic sites on the Euphrates using sexual dimorphism[END_REF][START_REF] Vigne | The early process of mammal domestication in the Near East: New evidence from the Pre-Neolithic and Pre-Pottery Neolithic in Cyprus[END_REF][START_REF] Zeder | Core questions in domestication research[END_REF].

The geometric morphometrics (GMM) revolution in the study of biological forms [START_REF] Zelditch | Geometric Morphometrics for Biologists: A Primer[END_REF] have provided an important breakthrough to track phenotypic responses associated with both the process of domestication and its geographical origins in the archaeological record [START_REF] Cucchi | Early Neolithic pig domestication at Jiahu, Henan Province, China: clues from molar shape analyses using geometric morphometric approaches[END_REF]Larson, Dobney & Cucchi 2011;[START_REF] Evin | A history of pig domestication: New ways of exploring a complex process[END_REF]. Mainly applied on suid teeth, GMM studies have contributed to reach new insights into pig domestication and dispersal in Island South East Asia [START_REF] Larson | Phylogeny and ancient DNA of Sus provides insights into neolithic expansion in Island Southeast Asia and Oceania[END_REF][START_REF] Cucchi | New insights into pig taxonomy, domestication and human dispersal in Island South East Asia: molar shape analysis of Sus remains from Niah Caves[END_REF], China [START_REF] Cucchi | Early Neolithic pig domestication at Jiahu, Henan Province, China: clues from molar shape analyses using geometric morphometric approaches[END_REF][START_REF] Cucchi | Social Complexification and Pig ( Sus scrofa ) Husbandry in Ancient China: A Combined Geometric Morphometric and Isotopic Approach[END_REF] and Europe [START_REF] Ottoni | Pig Domestication and Human-Mediated Dispersal in Western Eurasia Revealed through Ancient DNA and Geometric Morphometrics[END_REF][START_REF] Evin | Unravelling the complexity of domestication: a case study using morphometrics and ancient DNA analyses of archaeological pigs from Romania[END_REF][START_REF] Balasse | Wild game or farm animal? Tracking human-pig relationships in ancient times through stable isotope analysis[END_REF], as well as into the emergence of landrace breeds in Gaul [START_REF] Duval | Did Romanization impact Gallic pig morphology? New insights from molar geometric morphometrics[END_REF]2018) and Bronze-Age Switzerland [START_REF] Bopp-Ito | Phenotypic diversity in Bronze Age pigs from the Alpine and Central Plateau regions of Switzerland[END_REF]. The GMM approach of the enamel folding pattern of the occlusal surfaces (EFPOS hereafter) has been used to distinguish Bovidae from South Africa [START_REF] Brophy | Quantitative morphological analysis of bovid teeth and implications for paleoenvironmental reconstruction of Plovers Lake, Gauteng Province, South Africa[END_REF], explore diversity in Gallic-Roman domestic cattle [START_REF] Nuviala | Analyse des contours de la troisième molaire inférieure de boeuf : un moyen d'appréhender la variabilité morphologique du boeuf galloromain ?[END_REF] and assess taxonomic and phylogenetic signals in Equus mandibular teeth among extant wild, domestic and hybrid taxa [START_REF] Cucchi | Detecting taxonomic and phylogenetic signals in equid cheek teeth: towards new palaeontological and archaeological proxies[END_REF].

In this study we would like to further assess the biosystematic resolution of EFPOS variation in Bovini. The first objective is to test whether EFPOS can reach an accurate taxonomic identification at the genus and species levels to document the diversity in wild cattle before the emergence of agriculture. The second is to test whether genetic isolation can lead to observable intraspecific phenotypic divergence. And finally, we tested the extant of phylogenetic signal in EFPOS variation to assess its potential in tracking genetic relationship.

To provide a proof of concept useful across the wide geographic range where these questions apply, we have chosen to test these biosystematic signals from three molars (M 2 , M 3 and M3) commonly found in archaeological deposits and over two subtribes of Bovini: the Bubalina (Bubalus and Syncerus) and the Bovina (Bos and Bison) [START_REF] Hassanin | Combining multiple autosomal introns for studying shallow phylogeny and taxonomy of Laurasiatherian mammals: Application to the tribe Bovini (Cetartiodactyla, Bovidae)[END_REF], which covers most of the current diversity in wild and domestic cattle. We assessed which molar would perform best in identifying the genus and species of the Bovini tribe and populations within Bos p. taurus, taking into account the age and different parameters of the EFPOS quantification. We assessed the phylogenetic signal using a mtDNA genome phylogeny from modern specimens contrasted with the interspecific morphospace of each molar.

Materials and methods

Specimens analysed include three extant species of the subtribe Bubalina: the Asian water buffalo (Bubalus arnee), the lowland anoa (Bubalus depressicornis) and the African buffalo (Syncerus caffer); two extant species of Bison: the American bison (Bison bison), and the European bison (Bison bonasus); four wild species of the genus Bos: the yak (Bos mutus), the gaur (Bos gaurus), the banteng (Bos javanicus) and the extinct aurochs (Bos primigenius); and two domestic species: common cattle (Bos primigenius taurus) and zebu cattle (Bos primigenius indicus). The modern breeds of common cattle used in this study (Hollaender, Jersey, Shorthorn and Wilstermarsch) are from the Museum of Domesticated Animals "Julius Kühn" in Halle, which collected the specimens in the late 19 th and early 20 th century. These 167 specimens belong to the so-called ancient breeds or little developed country breeds. 168

Wilstermarsch and Hollaender belong to the Pied Lowland/Northwestern Lowland subunit, 169 whereas Shorthorn and Jersey belong to the British/Celtic/North European subunit (Felius et 170 al., 2011). Jersey and Hollaender are two milk breeds, while Wilstermarsch and Shorthorn are 171 dual-purpose/meat breeds. These different selective purposes are observable from their skull 172 morphotypes: milk breeds have longer and narrower skulls, whereas meat animals have 173 shorter and more compact ones [START_REF] Duerst | Grundlagen der Rinderzucht: eine Darstellung der wichtigsten für die Entwicklung der Leistungen und der Körperformen des Rindes ursächlichen, physiologischanatomischen, zoologisch-paläontologischen, entwicklungsmechanischen und kulturhistorischen Tatsachen und Lehren : für wissenschaftliches und praktisches Studium[END_REF][START_REF] Haring | Mast-und Schlachteigenschaften: und ihre Beziehungen zum Typ verschiedener Schweinerassen und deren Kreuzungen[END_REF]. We also included two 174 archaeological samples of domestic cattle (Bos p. taurus), one from the late Iron Age site of 175 Switzerland, and one from the Bronze Age site of 176 Savognin-Padnal (GR, Switzerland, c. 1950-1050 BC). 177 Due to a nomenclature problem (e.g. [START_REF] Gentry | The naming of wild animal species and their domestic derivatives[END_REF][START_REF] Roskov | Species 2000 & ITIS Catalogue of Life, 2018 Annual Checklist. Digital resource at[END_REF] Tab. 1: Sample size for the modern and archaeological Bovini used in this study. Maxillary molars: M 3 , M 2 ; mandibular molar: M3.

Our morphometric study focused on the maxillary second and third molars (M 2 , M 3 ) and the mandibular third molar (M3) (fig. 2a,2b). Both maxillary and mandibular third molars provide important advantages when dealing with archaeological material since they are easily distinguishable from the other teeth. Understanding the effect of attrition through ageing over the variation of the EFPOS was pivotal for this study. This effect has been proven null or negligible on the EFPOS of horses' cheek teeth [START_REF] Seetah | A geometric morphometric reevaluation of the use of dental form to explore differences in horse (Equus caballus) populations and its potential zooarchaeological application[END_REF], and CT analysis on the third maxillary and mandibular molars of African bovids [START_REF] Brophy | Quantitative morphological analysis of bovid teeth and implications for paleoenvironmental reconstruction of Plovers Lake, Gauteng Province, South Africa[END_REF]) also indicate a negligible effect. To investigate whether this observation is also valid for bovine teeth, we only considered teeth with full abrasion for this study. To obtain a relative age determination for the modern and archaeological specimens, we used the wear stage recording systems accessible from the IPAS website (https://duw.unibas.ch/de/ipna/forschung/archaeobiologie/archaeozoologie/methodik/). With this system, each tooth is given a numerical score which corresponds to its eruption and attrition stages, combining the age systems of [START_REF] Barone | Anatomie comparée des mammifères domestiques. Tome premier: Osteologie[END_REF][START_REF] Becker | Tierknochenfunde, zweiter Bericht: Mittleres und oberes Schichtpaket (MS und OS) der Cortaillod-Kultur. Die neotlithischen Ufersiedlungen von Twann[END_REF], Ducos 1968[START_REF] Habermehl | Die Altersbestimmung bei Haus-und Labortieren, 2[END_REF][START_REF] Grigson | Sex and age determination of some bones and teeth of domestic cattle: a 27 review of the literature[END_REF].

Geometric morphometric approach to Bovini molar's EFPOS

Following previous studies on equines [START_REF] Cucchi | Detecting taxonomic and phylogenetic signals in equid cheek teeth: towards new palaeontological and archaeological proxies[END_REF], we quantified the EFPOS of bovine molars with a 2D GMM approach combining landmarks of type 3 on the main enamel folds and semilandmarks on enamel curves between each landmark. For all three molars, we collected 7 landmarks (1-7) on the maximum curvature of the enamel folds at the junction between dentine and enamel, (fig. 2AB). The semilandmarks were then sampled equidistantly on the curves drawn between each landmark, at the junction between dentine and enamel (tab.

2). We performed a Generalised Procrustes Superimposition (GPS) on all the 2D Cartesian point coordinates of the EFPOS outline configurations to remove the effects of position, size and orientation. The semilandmarks were slided using the Bending Energy method (Bookestein 1996). The GPS produces Procrustes coordinates which can be used as EFPOS' shape parameter. GPS also computes the centroid size of the EFPOS as the square root of the sum of the squared distances from the landmarks and semilandmarks to the centroid of the enamel outline [START_REF] Bookstein | Landmark methods for forms without landmarks: morphometrics of group differences in outline shape[END_REF][START_REF] Zelditch | Geometric Morphometrics for Biologists: A Primer[END_REF] The points and curves digitization were computed from the images using tpsDig V2.17 and the GPS with tpsRelw v1.53 [START_REF] Rohlf | Morphometrics at SUNY Stony Brook[END_REF].

Assessing taxonomic signals in bovine molars EFPOS size and shape while accounting for age effect and allometry

Differences in centroid size variation of the EFPOS among bovine taxa were tested for the three molars using Analysis of Variance (ANOVA) with Bonferroni correction for pairwise comparisons (SI1). EFPOS size differences between bovine species, while accounting for size varying with age, were tested for each molar with a single factor Analysis of Covariance (ANCOVA). EFPOS shape differences between bovine species, while accounting for shape covarying with age and centroid size (allometry), were tested for each of the three molars with a single factor Multivariate Analysis of Covariance (MANCOVA). The MANCOVA assessed whether age and allometry effects differ significantly between the bovine species (interaction) and which of the three molars was the most affected by these two effects (R-squared).

We then assessed which of the three phenotypic markers (M 3 , M 2 , M3) would provide the greatest taxonomic resolution at the genus, species and population level, and which parameters of EFPOS quantification parameters (the shape, the size-free shape, the age-free shape or the form) would perform the best. The size-and age-free shape parameters are the residuals of the multivariate regression between the Procrustes coordinates as dependent variables and respectively the Log-centroid size and the age scores. The Procrustes from space were obtained by combining the Log-centroid size with the Procrustes coordinates [START_REF] Mitteroecker | Advances in Geometric Morphometrics[END_REF]Klingenberg 2015).

The taxonomic resolution (genus, species, and populations) for the three molars was assessed using Linear Discriminant Analysis (LDA). Its predicting approach calculated the percentage of cross validated correct classification (% CVC) of each specimen into its taxonomic group, following a 2 fold cross-validation over 10 000 iterations.

The genus factor comprised four groups: Bos, Bison, Bubalus and Syncerus. The species factor comprised the 11 specific groups of our dataset (tab. 1). The population factor was characterized by the 4 modern breeds and 2 archaeological populations (tab. 1). To avoid spurious results of the LDA induced by the high dimensionality of the shape parameters compared to the number of observation, we performed Principal Component Analysis (PCA) for each parameter and used only the PC scores explaining 99% of the total variance as variables in the LDA analysis. On average, 30 PC scores were enough to reach this threshold.

Assessing the phylogenetic signal in bovine molars EFPOS

To assess the phylogenetic signal of the three molars morphospace, we used the sequences of the complete mitochondrial genome of Bovini, extracted from the NCBI (National Center for Biotechnology Information) for the 11 following taxa (the accession numbers are indicated between brackets): domestic Asian buffalo -Bubalus arnee bubalis (JN632607), anoa - The sequences were aligned on SeaView version 4 [START_REF] Gouy | SeaView version 4: a multiplatform graphical user interface for sequence alignment and phylogenetic tree building[END_REF]. Our final alignment contained 16,284 characters. A Bayesian tree was constructed with MrBayes v3.2.1 [START_REF] Ronquist | MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space[END_REF] and the GTR+I+G model selected under jModelTest 2.1.7 [START_REF] Darriba | jModelTest 2: more models, new heuristics and parallel computing[END_REF] using the Akaike information criterion. The posterior probabilities were calculated using four independent Markov chains run for 10,000,000 Metropolis-coupled MCMC generations, with tree sampling every 1000 generations and a burn-in of 25%.

To estimate the strength of the phylogenetic signal relative to a Brownian motion model of evolution, we used the Kmult method [START_REF] Adams | A method for assessing phylogenetic least squares models for shape and other high-dimensional multivariate data[END_REF], a multivariate generalization of the Kstatistic [START_REF] Blomberg | TESTING FOR PHYLOGENETIC SIGNAL IN COMPARATIVE DATA: BEHAVIORAL TRAITS ARE MORE LABILE[END_REF]. Under a Brownian motion model of evolution the expected Kmult value is 1.0. Observed Kmult values were evaluated by comparing with the permutated K-value from randomized shape data relative to the tree topology of the molecular phylogeny.

To visualise the association between the molecular phylogeny and the dental shape space we created a phylomorphospace, which fitted the shape changes from the bovini taxa onto the tree topology using a squared change parsimony approach [START_REF] Klingenberg | Testing and quantifying phylogenetic signals and homoplasy in morphometric data[END_REF].

Morphometric statistics were performed with R (RCoreTeam, 2016), and the libraries MASS [START_REF] Venables | Modern applied statistics[END_REF] and geomorph [START_REF] Adams | geomorph: an R package for the collection and analysis of geometric morphometric shape data[END_REF].

Results

EFPOS size and shape differences among Bovini accounting for age and allometric effect

The ANCOVA found significant overall differences in the size variation of the EFPOS among the Bovini species (modern and archeological Bos p. taurus have been separated) for each tooth (Figure 3, tab. 3, See SI1 for taxa pairwise comparisons for each molar). All three molars displayed the same pattern of size differentiation among taxa, although the differentiations had a greater resolution for the maxillary molars. Bos p. primigenius is significantly bigger than all the other species, while Bubalus depressicornis is the smallest.

After Bos p. primigenius, the bisons' M 3 and M 2 were significantly larger than those of other genera. We also observed that Iron Age common cattle were significantly smaller than modern cattle breeds. In most cases there was an overlap in the size range of the remaining taxa. The size differences among Bovini species are reflected in the strong taxonomic signal in the size variation of the EFPOS which is above 50% (R-squared) for all three molars, while an age effect below 10% was only significant for the third maxillary and mandibular molars (tab.3), common among the species (tab. 3: Age*Taxa). MANCOVA found significant overall EFPOS shape differences among the Bovini species for the three molars (tab. 4: "Taxa"). Patterns of shape variation differ according to the molars (figure 4), but most of the variation was triggered by the differentiation between the Bubalina For both third molars, the age related attrition effect accounted for 10% and 12% of the EFPOS shape variation, while for the second maxillary molar the age effect reach 32% (tab. 4). It was reduced to around 18% by eliminating the youngest and oldest age groups. A common, and almost common, age effect among species for the third maxillary lower and mandibular molars for the M 2 (tab.4) was found. We then obtained age-free shape variables using the residuals of the multivariate regression. Tab. 4: Results of the single factor MANCOVA test for age effect on the shape variation of the M 3 , M 2 and M3, (Age effect) and the shape differences among Bovini species (Taxa) while accounting for shape covarying with age (Age*Taxa) , where"*" denotes interaction.

and
The MANCOVA found significant but very weak allometric effects in the M 3 and M3 shape variation, explaining less than 5% of the shape variation (tab. 5). Furthermore, this allometric effect did not interact with the taxonomic signal, which suggests a common allometric effect on molar shape that does not affect species differentiation. We obtained size-free shape variables using the residuals of the multivariate regression. Table 5: Results of single factor MANCOVA test for the allometric effect over shape variation of the M 3 , M 2 and M3 as well as shape differences between Bovini species (Taxa) while accounting for shape covarying with size (Log CS), where"*" denotes interaction.

Taxonomic signals in Bovini cheek teeth shape and form

The taxonomic resolution from the genus to the population levels differed among the phenotypic markers and their morphometric variables, according to the percentage of CVC (tab. 6). Maxillary molars performed better than the mandibular ones in separating the taxa to the genus and species level, which is not surprising since only half of the M3 could be quantified. The M 2 shape was the most efficient at discriminating the different genera in bovines with almost 90% correct classifications, closely followed by both third molars.

Discrimination of species was still high for the M 3 and M 2 shapes with more than 75% of correct classifications for the M 3 , but was enhanced in all three molars when using the form variables with M 2 , reaching nearly 85% CVC. For population discrimination the CVC dropped below 50%, with the M3 providing the best signal (48%). Finally, correcting for age or size effects did not significantly improve the discrimination of genus and species in bovines. Tab. 6: Percentage of correct classification after cross validation for M 3 , M 2 , and M3 at different taxonomic levels: genus, species and population and according to four phenotypic parameters: Shape, Size-free shape, Age-free shape and Form. *: analysis performed on extant populations of Bos p. taurus only.

Shape

CVA shows significant divergence (MANOVA: Pillai=1. 6655, F=2.9535, p<0.0001) between the ancestral wild shape of the Aurochs and the domestic phenotype of common and zebu cattle, between common and zebu cattle and between modern and archaeological common cattle, with zebu cattle showing the greatest divergence (Figure 5A). The CVC is too low to accurately separate the different populations of modern cattle based on their molar shape (Tab. 6). However, we found significant differentiation patterns (MANOVA: Pillai=2.423,

Phylogenetic signal in Bovini molar shape.

The Kmult statistics found a low but significant phylogenetic signal in M 3 (Kmult observed = 0.206; P-value = 0.009) and M3 (Kmult observed = 0.2781; P-value = 0.003) but not in the M 2 (Kmult observed = 0.17

Discussion

Identifying the diversity of the wild progenitors and the domestication process in Bos species

Our study has shown that both the EFPOS size and shape variations of the maxillary molars provided reliable taxonomic resolution at the genus and species levels. The resolution to species was not affected by the age effect, or only slightly for the M 2 shape, and could be increased when using the form parameter, bringing together the taxonomic resolution of both size and shape. These results are congruent with previous studies on South African Bovids [START_REF] Brophy | Quantitative morphological analysis of bovid teeth and implications for paleoenvironmental reconstruction of Plovers Lake, Gauteng Province, South Africa[END_REF]) and equids [START_REF] Cucchi | Detecting taxonomic and phylogenetic signals in equid cheek teeth: towards new palaeontological and archaeological proxies[END_REF] in rejecting the general idea that hypsodont teeth from ruminants are not good taxonomic characters because of the attrition blurring process. As expected the second molar was the most affected by attrition, leading to potential spurious identification. Therefore, we suggest to implementing taxonomic studies on the M 3 .

Our study also showed that both maxillary and mandibular third molar shape variations had a significant phylogenetic signal among Bovini species, finding phenotypic relationships congruent with the most recent phylogeny [START_REF] Hassanin | Systematic and evolution of Bovini[END_REF][START_REF] Hassanin | Systematics and phylogeny of cattle[END_REF].

Both molars captured the genetic relationship between aurochs (Bos p. primigenius) and the two lines of domestic cattle (Bos p. taurus / indicus), between the two species of bison (Bison bison and Bison bonasus) and the yak (Bos mutus), and between gaur (Bos gaurus) and banteng (Bos javanicus). However, the strength of this signal was weakened by the phenotypic relationship between the south-east Asian cattle (Bos gaurus and Bos javanicus)

and the Bubalina, as well as the one between Bos mutus and Bos p. indicus for the M 3 and the phenotypic proximity between Bos javanicus and Bos p. taurus for the M3. Since we could not make a distinction between the wild and domestic forms of yak (Bos mutus), gaur (Bos gaurus), and banteng (Bos javanicus) for this analysis, it is likely that this incongruence between genotype and phenotype is the consequence of genetic introgression among Bos species [START_REF] Hassanin | Systematics and phylogeny of cattle[END_REF][START_REF] Wu | Pervasive introgression facilitated domestication and adaptation in the Bos species complex[END_REF].

The high taxonomic resolution up to the species level, the lack of significant impact of attrition in this resolution and the strength of the phylogenetic signal recorded in the M 3 clearly suggest that this molar should be preferred for biosystematic studies of Bovini dental series from fossil and archaeological contexts. Since the wild ancestors of domestic cattle are either extinct (aurochs) or reduced to very small pockets of range distribution, this biosystematic power would prove extremely valuable to document the diversity and distribution of the wild progenitors of domestic forms, from the Pleistocene and early The biosystematic resolution of the M 3 could also prove useful to document the subsequent dispersal of domestic forms outside the distribution range of their wild ancestor.

For example, the possibility of identifying the two domestic cattle (Bos p. taurus and indicus)

and their relationships with their potential wild progenitor could help provide new clues on the development of herding in Africa, imported from the Near East during the 7 th millennium BC [START_REF] Linseele | Early Stock Keeping in Northeastern Africa: Near Eastern Influences and Local Developments[END_REF][START_REF] Lesur | The Advent of Herding in the Horn of Africa: New Data from Ethiopia, Djibouti and Somaliland[END_REF]. According to the latest understanding, the first domestic bovine in Africa are considered to be taurine cattle -zebu seem to occur later, during the 2 nd millennium BC -but we still do not know where it entered the continent, although genetic data suggest it happened through the Horn of Africa [START_REF] Hanotte | African Pastoralism: Genetic Imprints of Origins and Migrations[END_REF][START_REF] Gifford-Gonzalez | Domesticating Animals in Africa: Implications of Genetic and Archaeological Findings[END_REF].

The scarcity of evidence for zebu in archaeological contexts partly lies in the fact that only limited morphological criteria have been established to identify zebu from taurine specimens [START_REF] Grigson | The Craniology and Relationships of Four Species of Bos. 3. Basic Craniology: Bos taurus L. Sagittal Profiles and Other Non-Measurable Characters[END_REF][START_REF] Grigson | The Craniology and Relationships of Four Species of Bos. V. Bos indicus L[END_REF]. However, these criteria are not especially useful for archaeological material because the study was based on cranial parts that are seldom preserved in ancient contexts. For the moment, only a few studies have suggested the ancient presence of zebu in Africa, and only from the 1 st millennium BC [START_REF] Chaix | A Review of the History of Cattle in the Sudan throughout the Holocene[END_REF][START_REF] Chaix | The Fauna from the UNO/BU Excavations at Bieta Giyorgis (Aksum) in Tigray, Northern Ethiopia: Campaigns 1995-2003; Pre-Aksumite, 700-400 BC to Late 25 Aksumite, AD 800-1200[END_REF].

The GMM approach of EFPOS, therefore, provides a new method to address this issue.

However, the remaining difficulty lies in the fact that hybridization between the two species may have occurred very early in Africa [START_REF] Hanotte | African Pastoralism: Genetic Imprints of Origins and Migrations[END_REF][START_REF] Lesur | Et la gazelle devint chèvre. Pré-Histoires africaines d'hommes et d'animaux[END_REF]. The method now needs to integrate hybridised African specimens to check if we can also discriminate them.

Identifying the emergence of early cattle morphotype in archaeology

Identifying the emergence of early domestic cattle breeds due to human selective pressure, for traits such as greater milk yield, is in the realm of palaeogenomic studies [START_REF] Bradley | Genetics and the Origins of Domestic Cattle[END_REF], Irving-Pease et al. 2018) but not of morphometrics. Indeed, morphometric changes in skeleton bones cannot be directly related with a specific selective pressure [START_REF] Zeder | Core questions in domestication research[END_REF], such as meat or milk yield. However, documenting the origin and pace of morphometric changes over a large dental series from the 4 th millennium BC in the Near and Middle East and from the 3 rd millennium BC in Europe until Antiquity, would allow us to pinpoint the timing of the emergence of regional morphotypes due to neutral genetic shifts fostered by geographic isolation, as has already been performed for pigs (Duval et al. 2014[START_REF] Duval | The development of new husbandry and economic models in Gaul between the Iron Age and the Roman Period: New 26 insights from pig bones and teeth morphometrics[END_REF]. Once these key locations and periods are identified, targeted paleogenomic analysis could then be more efficiently implemented to explore selective traits within these regional morphotypes.

Our study found clear M 3 shape divergence between aurochs, Bronze and Iron Age Switzerland cattle and modern cattle breeds, suggesting that regional morphotypes are clearly attested in this region from the Bronze Age. Furthermore, these archaeological cattle morphotypes show greater phenotypic similarity towards modern cattle breeds than to their aurochsen progenitor, suggesting that great selective breeding in cattle was already in place from the 2 nd millennium BC in Europe. But was that selection directed towards adaptation to mountainous regions or towards greater meat or milk production? Only paleogenomic analyses would be able to provide a definite answer. The population signal in the molar shape variation of Bos p. taurus is too weak to make safe taxonomic identification. However, phenotypic proximities between the Bronze Age and Iron Age morphotypes with the Hollaender breed [START_REF] Felius | On the Breeds of Cattle -Historic and Current Classifications[END_REF], suggest that these Bronze and Iron Age morphoptypes may have been selected for milk production. Other potential explanations could be the genetic consequence of different regional connections between the two sites. Archaeological finds evidence that the Padnal people had connections with the north of Italy and with the north and east of Switzerland [START_REF] Rageth | Die wichtigsten Resultate der Ausgrabungen in der bronzezeitlichen Siedlung auf dem Padnal bei Savognin (Oberhalbstein GR)[END_REF]. The Basel-Gasfabrik people, on the other hand, had more connections with the west and south of France as well as with the Mediterranean [START_REF] Spichtig | Die spätlatènezeitliche Siedlung Basel-Gasfabrik[END_REF]. This hypothesis will have to be tested with a paleogenetic approach and a large comparative molar shape analysis of the Bronze and Iron Age domestic cattle sampled from Switzerland and Italy.

CONCLUSION

Our understanding of the origin and trajectories of wild cattle domestication have been mainly driven by the advance in genetic studies of modern and ancient Bos species, but zooarchaeologists still need powerful biometric markers to document the timing and mechanisms of cattle domestication in the old world and implement targeted biomolecular studies.

Our GMM approach of the taxonomic and phylogenetic signals encapsulated in the variation of molar morphology of wild and domestic species of the tribe Bovini found that M 3 shape variation could accurately separate the genera and species of both wild and domestic modern specimens of the tribe Bovini with high statistical accuracy, as well as their phylogenetic relationships. But recent genetic introgression between wild and domestic Bos species in eastern Asia have partly blurred this phylogenetic signal. The taxonomic resolution at a population level does not reach such a reliable probabilistic classification. However, significant differentiation patterns among aurochs and modern and ancient cattle breeds clearly suggest that genetic divergence indirectly triggered by intentional breeding selection can be deciphered from the archaeological record.

This proof of concept strongly suggest that GMM analysis of M 3 applied on Bos remains from Pleistocene and Holocene dental series and through Eurasia would be a useful pre requisite for targeted bio molecular studies for three main reasons: (1) to document the diversity and distribution of wild progenitors of domestic forms; (2) to identify where and when the first domestic phenotype emerged; and (3) to trace the emergence of the earliest regional morphotype and their trajectory towards modern breeds.

Figure 1 :

 1 Figure 1: Former range (Late Pleistocene) and current geographic distribution of wild species of Bovina after Hassanin 2015 with localisation of their potential domestication centres after Irving-Pease et al. 2018.

  Fig. 2: A: Landmarks (1-7) and semilandmarks (a-f) used for the maxillary M 2 and M 3 (here a M 3 of an archaeological Bos p. taurus from the Basel-Gasfabrik sample). B: Landmarks (1-7) and semilandmarks (a -f) used for the M3 (here a M3 of an archaeological Bos p. taurus from the Basel-Gasfabrik sample). Scale in mm.

  Bubalus depressicornis (EF536351), African buffalo -Syncerus caffer (EF536353), American bison -Bison bison (JN632601), European bison -Bison bonasus (JN632602), domestic yak -Bos mutus grunniens (EF494179), domestic cattle -Bos primigenius taurus (V00654), zebu -Bos primigenius indicus (AY126697), gaur -Bos gaurus (JN632604), banteng -Bos javanicus (JN632606) and aurochs -Bos primigenius primigenius (GU985279).

Fig. 3 :

 3 Fig. 3: Box plot of the M 3 , M 2 and M3 centroid sizes (log) from the tribe Bovini, including modern and archaeological domestic cattle (Bos p. taurus).

Figure 4 :

 4 Figure 4: Patterns of shape variation among the Bovini species mean shape for M 3 (A), M 2 (B) and M3 (C). Shape deformation along the PCs axes are displayed using the shape of specimens at the ends of the range of variability along PC1 and PC2.

F=1. 558

 558 Fig. 5: Patterns of M 3 shape differentiation within the first three canonical factors of the CVA (A) among wild extinct aurochs (Bos primigenius) and modern and archaeological domestic cattle (Bos p. taurus) and (B) among modern and archaeological common cattle.

  ; P-value = 0.18). Both third molars phylomorphospace captured the phylogenetic signal separating the Bovinae from the Bubalina and Bos from Bison species within the Bovina (fig. 6). The low phylogenetic signal in M 3 is driven by the phenotypic relationships of the cattle of south-east Asia -the banteng (Bos javanicus) and the gaur (Bos gaurus) within the Bubalina, and greater dental shape similarity between the yak (Bos mutus) and the humped cattle (Bos p. indicus) than with the American bison (Bison bison), both being incongruent with the mitochondrial phylogeny (fig. 6A). For the M3 (figure 6B), the low phylogenetic signal is due to greater phenotypic similarities of Bos p. primigenius and the common cattle (Bos p. taurus) with the American bison (Bison bison), than with the European one (Bison bonasus), as well as a greater similarity between the bateng (Bos javanicus) and the common cattle (Bos p. taurus) than with the gaur (Bos gaurus).

Fig. 6 :

 6 Fig. 6: Bovini MtDNA genome phylogeny and dental shape Phylomorphospace of M 3 (A) and

  Holocene archaeological records. Applied to all the available dental series of Bovini in the old world, this marker would help illustrate the natural geographic distribution of aurochs variants and test the relevance of the three different sub-species proposed by zooarchaeologists(Clutton-Brock 1989). In south-east Asia, these studies would help map the distribution of the four wild species of Bos and understand their geographic diversity and sympatry.The discrimination and the phylogenetic signal recorded in the M 3 between the wild and domestic European forms of Bovini (Bison bonasus, Bos p. primigenius and Bos p. taurus, Bos p. indicus) clearly suggest that GMM on M 3 could trace the phylogenetic histories of domestic cattle in the core area of the Fertile crescent. However, this approach alone could prove unable to identify the early process of domestication, when high gene flow between wild and domestic populations would prevent the expression of genetically induced phenotypic divergence(Irving-Pease et al. 2018[START_REF] Vigne | Early domestication and farming: what should we know or do for a better understanding?[END_REF]. Unless corralling and reproductive control were strong enough to foster the adaptive phenotypic change induced by the anthropic environment[START_REF] Zohary | The role of unconscious selection in the domestication of sheep and goats[END_REF]).

Savognin-Padnal (Bronze Age) 15

Acknowledgements

We are grateful to everyone who granted us access to museum collections and/or provided working space: L. Costeur (Naturhistorisches Museum Basel), M. Hellmund † and M. Stache (Geiseltalmuseum Halle), F. Mayer and C. Funk (Museum für Naturkunde Berlin) and K.

Murphy Gregersen (Zoologisk Museum Kopenhagen). For access to the archaeological samples we thank M. Bopp-Ito (IPAS Basel) and T. Reitmaier (Archaeological Service of the canton Graubünden), and G. Lassau and N. Spichtig (Archaeological Service of the canton Basel-Stadt). We are most grateful to the two anonymous reviewers who provided many insightful comments which helped us to greatly improve the quality of this manuscript. Many thanks to Jill Cucchi for copy-editing.

Authors contributions:

TC and BS designed the study, BS collected the morphometric and Age data, TC and BS analyzed and interpreted the data, TC and BS wrote the paper in collaboration with RS, AH, JL and JS.

Funding:

Barbara Stopp has been supported by the Swiss National Science Foundation (SNSF: PP00P1_123336) and the Institute of Prehistory and Archaeological Science (IPAS) at the University of Basel. This study did not receive any other specific grants from funding agencies in the public, commercial or not-for-profit sectors.


T R Gimadiev

O Klimchuk

R I Nugmanov

T I Madzhidov

A Varnek

Sydnone-alkyne cycloaddition: which factors are responsible for reaction rate ?

Keywords: biorthogonal reactions, sydnones, DFT, reaction rate, structurereactivity models

In this work we report extensive DFT study of sydnone-alkyne cyclization which included investigation of the reaction mechanism, analysis of different factors affecting sydnone and alkyne reactivity as well as attempt to reproduce quantitatively experimental activation free energy. The calculations were performed for a set of 18 sydnone-alkyne reactions with a help of a semiautomatized workflow involving reagent preparation and generation of starting structures for a plausible transition state. Reconstructed reaction path supported two-step mechanism: cycloaddition followed by retro-Diels-Alder reaction. Since the latter had a tiny barrier, the cycloaddition step was predicted to be the ratelimiting. For the ensemble of reactions, calculations reproduce activation free energies extracted from experimental reaction rates (k) with the accuracy of 2 kcal/mol. Accounting for solvation effects didn't change the overall trend of activation free energies as a function of substituents. A series of statistical model linking logk and sydnones structure was built using Support Vector Regression and Multiple Linear Regression machine-learning methods coupled with different types of molecular descriptors; none of them demonstrated a good performance at crossvalidation stage. Detailed analysis of different factors affecting reaction rate variation as a function of substituents revealed particular role of the charge on C3 atom in the sydnone moiety as well as of the size of the substituent at C3.

Introduction

The ability to selectively form and break chemical bonds in chemically complex and uncontrollable biological media is a long-standing goal of chemists interested in modifying biological materials. Biorthogonal chemical reactions [START_REF] Sletten | From Mechanism to Mouse: A Tale of Two Bioorthogonal Reactions[END_REF],

which don't interfere with biological processes, address precisely this challenge and therefore are of major importance in the fields of chemical biology and biochemistry. To fulfill the requirements of bioorthogonality, reaction partners must be stable and inert towards the plethora of chemical functionalities found in living systems while reacting selectively, efficiently and rapidly under physiological conditions with no or innocuous by-products [START_REF] Lim | Bioorthogonal chemistry: recent progress and future directions[END_REF][START_REF] Gong | Recent advances in bioorthogonal reactions for site-specific protein labeling and engineering[END_REF]. The coppercatalyzed 1,3-dipolar cycloaddition reactions, e.g., azide-alkyne cycloaddition, are the most extensively studied biorthogonal transformations. However, the Cu salts used as catalyst are very toxic for living organisms [START_REF] Kennedy | Cellular Consequences of Copper Complexes Used To Catalyze Bioorthogonal Click Reactions[END_REF]. This was a reason, to explore new reactions based on strain promoted cycloadditions [START_REF] Mbua | Strain-Promoted Alkyne-Azide Cycloadditions (SPAAC) Reveal New Features of Glycoconjugate Biosynthesis[END_REF] which proceed without any catalyst but, unfortunately, have rather low rate [START_REF] Sletten | Bioorthogonal Chemistry: Fishing for Selectivity in a Sea of Functionality[END_REF].

Recently sydnones were discovered to be a promising class of biorthogonal 1,3-diene which can participate in cyclization reaction with cyclic alkynes resulting in pyrazole formation [START_REF] Browne | A Sydnone Cycloaddition Route to Pyrazole Boronic Esters[END_REF]. This reaction proceeds in two step: cyclization with formation of reactive intermediate and retro-Diels-Alder reaction leading to CO 2 cleavage (Figure 1). Thermal cycloaddition of sydnones with alkynes require quite harsh conditions and proceeds with low regioselectivity [START_REF] Browne | A Sydnone Cycloaddition Route to Pyrazole Boronic Esters[END_REF] while copper catalysed reaction proceed selectivity in mild conditions with high yield [START_REF] Kolodych | Discovery of Chemoselective and Biocompatible Reactions Using a High-Throughput Immunoassay Screening[END_REF][START_REF] Specklin | One-Pot Synthesis of 1,4-Disubstituted Pyrazoles from Arylglycines via Copper-Catalyzed Sydnone-Alkyne Cycloaddition Reaction[END_REF][START_REF] Decuypere | Copper(I)-Catalyzed Cycloaddition of 4-Bromosydnones and Alkynes for the Regioselective Synthesis of 1,4,5-Trisubstituted Pyrazoles[END_REF]. Strained alkynes readily react with sydnones and iminosydnones without catalyst [START_REF] Wallace | Strain-promoted sydnone bicyclo-[6.1.0]-nonyne cycloaddition[END_REF][START_REF] Narayanam | Discovery of new mutually orthogonal bioorthogonal cycloaddition pairs through computational screening[END_REF][START_REF] Bernard | Bioorthogonal Click and Release Reaction of Iminosydnones with Cycloalkynes[END_REF].

Substituents in sydnone ring [START_REF] Plougastel | 4-Halogeno-sydnones for fast strain promoted cycloaddition with bicyclo-[6.1.0]-nonyne[END_REF][START_REF] Liu | Ultrafast Click Chemistry with Fluorosydnones[END_REF] as well as alkyne structure [START_REF] Narayanam | Discovery of new mutually orthogonal bioorthogonal cycloaddition pairs through computational screening[END_REF] greatly impact the reaction rate switching it from moderate to ultra-fast. Recently, detailed stopped-flow kinetics study of reaction of some sydnones supported two-stage mechanism [START_REF] Liu | Ultrafast Click Chemistry with Fluorosydnones[END_REF].

So far, theoretical studies were mostly focused on structure of alkynes or alkenes possessing the greatest reactivity in 1,3-cycloaddition. Thus, Narayanam et al [START_REF] Narayanam | Discovery of new mutually orthogonal bioorthogonal cycloaddition pairs through computational screening[END_REF] performed a DFT study of N-phenyl sydnone and diaryl-1,2,4,5-tetrazine reaction with eight different alkenes and alkynes. It has been demonstrated that sydnone interactions with the most strained alkynes proceed with large reaction rate. The importance of steric strain of alkene or alkyne moieties to accelerate kinetics of 1,3-cycloaddition with organic azides or tetrazines was also demonstrated in [START_REF] Lopez | Alkene Distortion Energies and Torsional Effects Control Reactivities, and Stereoselectivities of Azide Cycloadditions to Norbornene and Substituted Norbornenes[END_REF][START_REF] Liu | Theoretical Elucidation of the Origins of Substituent and Strain Effects on the Rates of Diels-Alder Reactions of 1,2,4,5-Tetrazines[END_REF]. It looks like strained alkenes and azides are well preorganized toward the Diels-Alder transition structure and, therefore, interact easier.

Tao et al. [START_REF] Tao | Origins of halogen effects in bioorthogonal sydnone cycloadditions[END_REF] performed detailed analysis of cycloaddition rate enhancement for halogen-substituted sydnones. They revealed that distortion energies of sydnone moiety well correlated with activation energies. In such a way, high reactivity of halosydnones was explained by their low distortion energies.

Analysis by Tao et al. [START_REF] Tao | Origins of halogen effects in bioorthogonal sydnone cycloadditions[END_REF] was limited to halogen-substituted sydnones and can't be directly used to predict reaction rate of cycloaddition involving sydnones with any substituents. Moreover, Tao et al. [START_REF] Tao | Origins of halogen effects in bioorthogonal sydnone cycloadditions[END_REF] considered only the first transition state of the reaction between sydnones and bicyclo[6.1.0]nonyne (BCN) but not retro-Diels-Alder reaction.

In this paper, we describe DFT-based full pathway of sydnone-alkyne reactions considering different sydnones and two alkynes -BCN and tetramethylthiacycloheptyne (TMTH). Our aim was to develop a fast semiautomatized screening workflow aiming to select reactants with desired reactivity.

A linear correlation between predicted and experimental activation free energies for 18 sydnone-alkyne reactions was observed which shows that gas-phase DFT calculations can be used as a screening tool for reaction design. We also attempted to build structure-reactivity models linking reaction rate with sydnones structure which provides an insight of different factors responsible for sydnones reactivity.

Computational details

Experimental reaction rate data

The data on cycloaddition reaction rate were taken from publications of F.

Taran's group [START_REF] Kolodych | Discovery of Chemoselective and Biocompatible Reactions Using a High-Throughput Immunoassay Screening[END_REF][START_REF] Bernard | Bioorthogonal Click and Release Reaction of Iminosydnones with Cycloalkynes[END_REF]. The collected dataset contains rate constants at ambient temperatures of 15 reactions of sydnones with BCN and 3 reactions of sydnones with TMTH, see Table 1 and Figure 1. In the series of reactions involving BCN, the fastest is reaction 15 where sydnone bears F substituent at C 3 position. As it follows from Table 1, the same sydnone interacts with TMTH much faster than with BCN (see reactions 15 and 16). Notice that second-order rate constants given in Table 1 correspond to cycloaddition step only because they were obtained by HPLC monitoring of sydnone conversion [START_REF] Plougastel | 4-Halogeno-sydnones for fast strain promoted cycloaddition with bicyclo-[6.1.0]-nonyne[END_REF].

Figure 1. The two steps reactions considered in this study: (3+2) cyclisation followed by retro-Diels-Alder reaction.

Quantum chemical calculations

Density Functional Theory (DFT) [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF][START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF] calculations were performed with the Priroda 11 program [START_REF] Laikov | PRIRODA-04: a quantum-chemical program suite. New possibilities in the study of molecular systems with the application of parallel computing[END_REF]22] using PBE exchange and correlation functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] and built-in triple-zeta split valence basis set 3z which is equivalent to Schäfer's TZVP basis set [START_REF] Schäfer | Fully optimized contracted Gaussian basis sets of triple zeta valence quality for atoms Li to Kr[END_REF]. Relativistic effects were neglected. Priroda11 is one of the fastest DFT software due to efficient evaluation of density functional exchange-correlation terms based on the electron density expansion [START_REF] Laikov | Fast evaluation of density functional exchange-correlation terms using the expansion of the electron density in auxiliary basis sets[END_REF].

Geometry optimization for reagents, intermediates and products as well as saddle point optimization were performed using built-in quasi-Newton-Raphson method and BFGS hessian update scheme. Scanning procedure was used to localize reasonable structural guess for transition state. This was achieved in constrained local optimization along the N1…C alkyne coordinate which varied from 2.1 to 3.1 Å with step 0.1 Å. In order to detect different kind of species (stable molecules, intermediates, transition state) geometry optimization was followed by the frequencies calculations. Thus, no imaginary frequencies must be found for stable species (reagents, products and intermediates) and only one large imaginary frequency must be detected for a transition state. The correctness of transition states was also checked using the intrinsic reaction coordinate (IRC) following approach.

Full reaction pathways were calculated for reaction 2 (also considered in early theoretical calculations [START_REF] Narayanam | Discovery of new mutually orthogonal bioorthogonal cycloaddition pairs through computational screening[END_REF]) as well as reactions 15 and 16 for which detailed kinetic experiments were performed by Liu et al [START_REF] Liu | Ultrafast Click Chemistry with Fluorosydnones[END_REF].

Since continuum solvation is not available in PRIRODA, for some molecules and transition states, solvation free energy in water was calculated using the IEF-PCM model [START_REF] Tomasi | The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level[END_REF][START_REF] Tomasi | Quantum Mechanical Continuum Solvation Models[END_REF] with SMD parameters for non-electrostatic terms [START_REF] Cramer | An SCF Solvation Model for the Hydrophobic Effect and Absolute Free Energies of Aqueous Solvation[END_REF] implemented in the Gaussian 09 program [START_REF] Frisch | Gaussian 09 Revision D.01[END_REF]. In this case, geometry optimization and hessian were calculated with Gaussian 09 [START_REF] Frisch | Gaussian 09 Revision D.01[END_REF] using PBEPBE/6-311++G(d,p) method.

Workflow for transition state detection

All structures were drawn with the Marvin View program [START_REF] Chemaxon | [END_REF] following the same numeration of atoms in sydnone ring. ChemAxon Calculator Plugin cxcalc [31] was used to generate up to 500 conformers with diversity limit 0.1 Å followed by geometry optimization with the MMFF94 force field [START_REF] Halgren | Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94[END_REF][START_REF] Halgren | Merck molecular force field. II. MMFF94 van der Waals and electrostatic parameters for intermolecular interactions[END_REF][START_REF] Halgren | Merck molecular force field. III. Molecular geometries and vibrational frequencies for MMFF94[END_REF]. Then, geometry of each structure was optimized using the PM6 semi-empiric method [START_REF] Stewart | Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements[END_REF] implemented in the MOPAC 11 program [36]. Duplicated conformers (with RMSD <0.1 Å) were removed using an in-house script which aligns molecules according to principal components of inertia. Remaining conformers were ranked according to their energies, and then 100 most stable conformers were retained for further calculations. All structures issued from PM6 optimization were reoptimized with DFT PBE/3z. One lowest energy conformer was selected.

Four possible orientations of reagents were considered (Figure 2). They lead to formation of two pairs of enantiomeric 1 st transition states (TS1). Further, at the retro Diels-Alder stage, a racemic mixture of two enantiomers was formed.

Four possible orientations of reagents were generated using in-house Python script in a way that distance between reaction center atoms (C3…C10 and N1…C11 or C3…C11 and N1…C10, depending on orientation) was 3.1 Å, i.e.

much greater than those in TS1 (Figure 2). A reasonable guess for TS1 was found using scan procedure implemented in PRIRODA11 in which the distance between reaction center atoms varied from 3.1 Å to 2.1 Å with step size 0.1 Å. For a greatest energy structure, the hessian was calculated followed by saddle point optimization.

Theoretically, pair of enantiomeric structures of TS1, intermediate and TS2

should have the same energies. In practice, however, due to possibility to generate different conformations for the stereoisomers, the energies of TS1 resulted from 4 different orientations of reagents slightly differ. For further analysis the lowest energy transition state among the four structures was selected. Hessian for the selected structure was calculated and presence of one large imaginary frequency was checked. The correctness of transition state TS1 was checked using Intrinsic Reactional Coordinate following (IRC) procedure.

The second transition state TS2 was found after IRC downhill shift towards the products. Resulting intermediate structure was used as starting structure for scanning reaction coordinate (N1-O5 bond length in this case). We found that small step of 0.01 Å was required because TS2 looks like a small "hubble" on the potential energy surface followed by steep energy drop.

Detailed description of the workflow and corresponding Python scripts are given in Appendix 1 in Supporting Information. Formula (1) was used for calculation of activation free energy on the basis of quantum chemical calculations:

∆G ∆G ∆G ∆G (1)
Calculation of activation free energy based on experimentally measured reaction rate was performed using formula (3) derived from Eyring equation of the Transition state theory (2) :

∆ !

(2)

∆" #$ln * ) (3)
where k -rate constant; κ -transmission coefficient (here, κ=1); T -Temperature in Kelvin; K B -Boltzmann constant; ∆∆G ≠ -Difference of free energies of reagents and TS at 298.15 K; h -Planck constant; R -gas constant.

Conceptual DFT indices calculation

Several Conceptual DFT [START_REF] Geerlings | Conceptual Density Functional Theory[END_REF] indices were used in order to reveal structural factors responsible for reaction rate variation. Electrophilicity index, ω, measures the stabilization energy when system acquires an additional electronic charge [START_REF] Parr | Electrophilicity Index[END_REF]:

* + 2-
where μ is electronic chemical potential, and η -chemical hardness. The latter could be expressed in terms of HOMO (ε HOMO ) and LUMO (ε LUMO ) energies as: Charge transfer from alkyne to sydnone in the transition state was calculated as absolute value of the entire electric charge of one of the species.

+ . /

QSPR modeling

Quantitative Structure-Property Relationships (QSPR) for logk were established using Support Vector Regression (SVR) and Multiple Linear Regression (MLR) machine-learning methods as well as ISIDA fragment descriptors [START_REF] Varnek | Substructural fragments: an universal language to encode reactions, molecular and supramolecular structures[END_REF] or quantum mechanics descriptors generated for sydnones.

SVR models were built and validated using the ε-SVR algorithm implemented in the libSVM package [START_REF] Chang | LIBSVM: A Library for Support Vector Machines[END_REF] using ISIDA fragment descriptors. The modeling was performed using the evolutionary SVR optimizer [START_REF] Horvath | An Evolutionary Optimizer of libsvm Models[END_REF], which can be used to perform both descriptor space selection and optimization of the operational parameters (epsilon, kernel type, cost, gamma) of the SVR method.

The MLR models were built using CODESSA software [START_REF] Karelson | QSPR and QSAR Models Derived Using Large Molecular Descriptor Spaces. A Review of CODESSA Applications[END_REF]. Three different sets of descriptors were explored: (i) quantum-chemical CODESSA-PRO descriptors [START_REF] Karelson | QSPR and QSAR Models Derived Using Large Molecular Descriptor Spaces. A Review of CODESSA Applications[END_REF], (ii) descriptors based on Bader's quantum theory of Atoms in Molecules [START_REF] Bader | Atoms in Molecules: A Quantum Theory[END_REF] calculated with CODESSA-3 [START_REF] Katritzky | CODESSA 3 program[END_REF] program, and, (iii) conceptual DFT [START_REF] Geerlings | Conceptual Density Functional Theory[END_REF] indices calculated for key atoms of sydnone with Priroda 11 [START_REF] Laikov | PRIRODA-04: a quantum-chemical program suite. New possibilities in the study of molecular systems with the application of parallel computing[END_REF]22].

Predictive performance of the models has been estimated by root mean squared error and squared determination coefficient (R 2 ) calculated in 10-fold cross-validation. Notice that predictive performance is acceptable only for the models with R 2 >0.5.

Results and Discussion

Reaction pathway investigation

The full reaction path was explored for the reaction 2 and 15 (see Table 1)

sharing the same alkyne (BCN) and sydnones differed by the substituent X (X = H and F, respectively). The procedure included scanning along the reaction Figure 3 shows that the limiting step of reaction 15 is the first transition state (TS1), free energy of the intermediate is 21 kcal/mol lower than that of reagents.

The intermediate looks rather unstable since it is separated from products by transition state (TS2) with a tiny barrier (some 1 kcal/mol). Decomposition of the intermediate is exergonic by some 85 kcal/mol. The pathway of reaction 2 looks very similar to that of 15. Much smaller, compared to 15, rate constant of 2 (Table 1) can be explained by ~2 kcal/mol higher activation free energy which is in line with the results by Tao et al [START_REF] Tao | Origins of halogen effects in bioorthogonal sydnone cycloadditions[END_REF]. Quite small TS2 activation free energies decrease from reaction 2 to 15 together with TS1 energies. Unlike reactions 2 and 15, reaction 16 between fluorosydnone and TMTH proceeds as synchronous one-step addition /elimination leading immediately to the products. Similar observations concerning very small second barrier in reactions between sydnone and BCN were also reported by Narayanam et al [START_REF] Narayanam | Discovery of new mutually orthogonal bioorthogonal cycloaddition pairs through computational screening[END_REF]. These observations, however, contradict with experimental results by Liu et al [START_REF] Liu | Ultrafast Click Chemistry with Fluorosydnones[END_REF] who demonstrated that the second step was rate limiting and quite slow compared to the first one. A contradiction between experiment and computations could be related to cage effect of solvent or some specific interactions TS2 with solvent. On the other hand, solvation free energies of reactants and TS1 estimated with the IEF-PCM model [START_REF] Tomasi | The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level[END_REF][START_REF] Tomasi | Quantum Mechanical Continuum Solvation Models[END_REF] were found very similar for all studied reactions. Thus, the first-step reaction barrier was found insensitive to solvent effects.

Activation free energies: prediction vs experiment

The question arises: whether rate constants estimated from activation free energy correlate with the experimentally measured values? Activation free energies, ΔG ≠ calc , were assessed for all 18 studied reactions according to eq. ( 1).

For the comparison purposes, "experimental" free energies, ΔG ≠ exp , were calculated from the rate constants according to eq. ( 3). Figure 4 demonstrates that calculated free energies well correlated with the "experimental" ones.

QSPR models

Quantum chemical assessment of reaction rate is rather reliable but, computational resources consuming task. In this case, QSPR models for logk of could become an efficient solution allowing to accelerate the screening of potentially interesting reactants. A large panel of models was built using Multiple Linear Regression and Support Vector Regression [START_REF] Drucker | Support vector regression machines[END_REF] machine learning methods and different types of molecular descriptors for sydnones including conceptual DFT indices (see Section 2.6 for details). Despite of our efforts, no models with acceptable predictive performance in cross-validation were obtained. This can be explained by small number of considered cycloaddition reactions and relatively large structural diversity of sydnones listed in Table 1.

Analysis of structural and electronic factors affecting reaction rate

Conceptual Density Functional Theory (DFT) indices [START_REF] Geerlings | Conceptual Density Functional Theory[END_REF] are often used to analyze cycloaddition reactions [START_REF] Domingo | Quantitative characterization of the global electrophilicity power of common diene/dienophile pairs in Diels-Alder reactions[END_REF][START_REF] Domingo | Understanding the mechanism of polar Diels-Alder reactions[END_REF]. To reveal structural factors affecting the rate constant we looked for linear correlations within restricted series of reactions involving the same alkyne (BCN). However, neither sydnone electrophilicity index [START_REF] Parr | Electrophilicity Index[END_REF] nor chemical potential (calculated as half-sum of HOMO and LUMO energies) nor chemical hardness (calculated as HOMO -LUMO gap) correlate with logk.

Analysis of orbital symmetries and energies shows that electrons are transferred from HOMO of alkyne to LUMO of sydnone (Figure 5) in agreement with earlier publication by Liang et al [START_REF] Liang | Control and Design of Mutual Orthogonality in Bioorthogonal Cycloadditions[END_REF] studied azide and tetrazine cycloaddition to alkynes. No valid correlation of LUMO of sydnones with reaction rate was found. In contrast to observations of Domingo et al [START_REF] Domingo | Understanding the mechanism of polar Diels-Alder reactions[END_REF] who studied Diels-Alder reactions between cyclopentadiene and some ethylene derivatives, no correlation of logk with the charge transfer in the transition state was found. So even within this congeneric series of reactions there is no evidence that orbital interactions represent the major factor affecting the reaction rate variations. Different descriptors characterizing electrostatic interactions (Fukui nucleophilicity, electrophilicity and radical attack susceptibility) using Conceptual DFT [START_REF] Geerlings | Conceptual Density Functional Theory[END_REF] were calculated. However, none of them correlated with logk values for all 15 considered sydnones -BCN reactions. One may suggest that there exists a complex interplay of structural and electronic factors which cannot be explained by simple linear correlations.

In order to reveal factors responsible for reaction rate, we considered two small congeneric series of reactions. The first one assembled reactions of BCN with sydnones having different substituents X and the same R = p-Me C 6 H 4 (2, 10-12, 15), whereas the second one included three reactions of TMTH with fluorosydnones (X=F) having different substituents R [START_REF] Lopez | Alkene Distortion Energies and Torsional Effects Control Reactivities, and Stereoselectivities of Azide Cycloadditions to Norbornene and Substituted Norbornenes[END_REF][START_REF] Liu | Theoretical Elucidation of the Origins of Substituent and Strain Effects on the Rates of Diels-Alder Reactions of 1,2,4,5-Tetrazines[END_REF][START_REF] Tao | Origins of halogen effects in bioorthogonal sydnone cycloadditions[END_REF]. Calculations for the first series show significant variation of charges on С3 atom, which change from negative (X = H) to positive value (X = F), whereas the charge on N2 atom slightly varies as a function of substituents. In this series, logk well correlate with the C3 charges (Figure 6a), r=0.98 and RMSE = 0.13 logk, p value 4•10 -4 . In line with observations reported by Tao et al. [START_REF] Tao | Origins of halogen effects in bioorthogonal sydnone cycloadditions[END_REF] , no correlation with LUMO energy was observed (Figure 6b).

Variation of aryl substituents at N3 atom in reaction of fluorosydnones with TMTH [START_REF] Lopez | Alkene Distortion Energies and Torsional Effects Control Reactivities, and Stereoselectivities of Azide Cycloadditions to Norbornene and Substituted Norbornenes[END_REF][START_REF] Liu | Theoretical Elucidation of the Origins of Substituent and Strain Effects on the Rates of Diels-Alder Reactions of 1,2,4,5-Tetrazines[END_REF][START_REF] Tao | Origins of halogen effects in bioorthogonal sydnone cycloadditions[END_REF], doesn't lead to significant variations of charge at C3 atom (Figure 6a). The LUMO energy in this series decreases from -3.57 to -3.05 eV which could partially explain considerable raise of the rate constant (Figure 6b). 1. Reactions with TMTH and BCN are shown in red and blue, respectively. The reaction numbers correspond to those in Table 1.

Notice that larger steric strain of alkyne in TMTH compared to BCN affects both charges on triple bond carbon atoms (-0.663 and -0.652 in BCN and -0.636 for both atoms in TMTH) and HOMO energies (-5.88 eV in BCN and -5.15 eV in TMTH, Figure 6). The latter favors larger reactivity of TMTH. This effect was mentioned in previous works on optimization of alkyne reactivity [START_REF] Liu | Bioorthogonal Cycloadditions: Computational Analysis with the Distortion/Interaction Model and Predictions of Reactivities[END_REF].

DFT calculations of the transition state of the cycloaddition of sydnones with BCN reveal an increase of the distance between alkyne C atom and C3 atom (C alkyne -C3) with the size of the substituent at C4 atom: 2.283 Å, 2.345 Å and 2.348 Å for X = H, F and I, respectively (Figure 7). This could be explained by steric Notice that steric and electronic effects for halogens vary in the same direction.

The sum of valence angles around C3 atoms can be considered as an implicit measure of its planarity: more its value differs from 360°, less planar (and more strained) the sydnone is. Notice that these observations disagree with conclusions drawn by Liang et al [START_REF] Liang | Control and Design of Mutual Orthogonality in Bioorthogonal Cycloadditions[END_REF] that electronically more reactive electrophiles should be sterically more encumbered. Indeed, iodo-sydnone in 12 is sterically more encumbered than fluoro-sydnone in 15. However, the former is less reactive due to weaker electron acceptor ability, on one hand, and steric repulsion with alkyne hydrogens, on the other hand. Replacement of F with more bulky electron-withdrawing substituent leads to reduction of the reaction rate. Thus reaction 9 involving sydnone with X=CF 3 is drastically slower than reaction 10 with X=Cl, despite the fact that atomic charges on C3 in these molecules are similar (0.009 and 0.0042 a.u., respectively).

In the context of rational design of highly reactive substrates, our calculations lead to the following observations:

• Strong electron acceptor X at C3 atom increases its positive charge and, thus, improves its affinity to alkyne in cycloaddition. However, this substituent should not be bulky because it may cause steric hindrance in the transition state. In this context, fluorine substituent at C4 represents an optimal choice.

• Strong electron acceptor at N3 atom weakly affects charge distribution in sydnone moiety but decreases its LUMO energy. This, in turn, favors cycloaddition, especially in reactions with TMTH.

• Compared to BCN, more sterically strained TMTH has higher HOMO energy which explains its better reactivity.

Conclusions

Extensive DFT calculations on a set of sydnone-alkyne reactions confirmed previously suggested two-step mechanism: cycloaddition followed by retro-Diels-Alder reaction. Since the latter had a tiny barrier, the cycloaddition step was predicted to be the rate-limiting. For the ensemble of 18 reactions, calculations reproduce activation free energies extracted from experimental reaction rates (k)

with the accuracy of 2 kcal/mol. Accounting for solvation effects didn't change the overall trend of activation free energies as a function of substituents. As an alternative, a series of statistical model linking logk and sydnones structure was built using Support Vector Regression and Multiple Linear Regression machinelearning methods coupled with different types of molecular descriptors; none of them demonstrated a good performance at cross-validation stage.

Detailed analysis of different factors affecting reaction rate variation as a function of substituents revealed particular role of the charge on C3 atom n the sydnone moiety (Figure 1) as well as of the size of the substituent at C3. Large positive charge on C3 enhances reaction rate due to transition state energy lowering, whereas bulky electron withdrawing substituent X at C3 atom leads to steric hindrance with alkyne hydrogens upon the transition state formation. Thus, only sterically unencumbered substituents X allow gaining reaction speed.

Analysis of orbital symmetries and energies shows that in the course of reaction, electrons are transferred from HOMO of alkyne to LUMO of sydnone.

The lower LUMO energy, the larger is reaction rate. Since LUMO energy is affected by π-electron withdrawing ability of substituent R, the usage of stronger electron acceptors will favor reaction.

Figure 2 .

 2 Figure 2. Four possible pre-reaction complexes sydnone and alkyne (top) and related products (bottom)

  coordinate, optimization of the transition state, motion along the intrinsic reaction coordinate which links transition state with reagents and products. At the initial stage, the first transition state (TS1) was localized. An intermediate was found in course of intrinsic reaction coordinate downhill movement with following local optimization. The intermediate's geometry was used to identify the second transition state (TS2) corresponding to retro-Diels-Alder reaction. In this case, an elongation of the N1-O5 bond was used as reaction coordinate (Figure2). Since TS2 structural was close to intermediate, a small step size of 0.01 Å was used in the scan procedure. Geometry of TS2 was found by saddle point optimization of the highest energy structure found in the scan procedure. IRC procedure supported that TS2 belong to reaction minimal energy pathway. Descent from TS2 resulted in either the intermediate or the product molecules. Energetic profile of reaction (Figure3) shows that the N1-O5 bond length smoothly elongates along the reaction path. Thus, for reaction 15 this bond length is equal to 1.379 Å in reagents, 1.414 Å in TS1, 1.546 Å in the intermediate and 1.721 Å in TS2.

Figure 3 .

 3 Figure 3. Reaction pathways for reactions 2 (dashed red line), 15 (solid black line) and 16 (dashed blue line). Relative free energies at 298K of molecules with respect to reagents are shown. Structure of reagents, transition states and intermediate for reaction 15 is shown. Substituent R of sydnone and some atoms of alkyne are omitted for the sake of clarity. Bond orders correspond to molecule representation on Figure 1.

Figure 5 .

 5 Figure 5. HOMO and LUMO energies obtained in DFT PBE/3z (TZVP) calculations for reactions 2, 15, 16. The energies are given in eV.

Figure 6 .

 6 Figure 6. Variation of logk in congeneric series of halo-sydnones as a function of (a) Hirshfield charge of C3 atom of sydnone and (b) LUMO of sydnone. The numbers correspond to reactions in Table1. Reactions with TMTH and BCN are shown in red and blue, respectively. The reaction numbers correspond to those in Table1.

  repulsion between substituent X and aliphatic CH 2 group next to triple bond in the alkyne moiety in the transition state. The shortest contact H…I (3.486 Å) in transition state TS1 is almost equal to the sum of the Pauling vdW radius of H and I (3.4 Å). While fluorine atom is almost located in the plane of sydnone ring in TS1, iodine atom deviates from it avoiding steric clashes, Figure 7. It can also be noticed that distances C3-C10 and N1-C11 in case of fluoro-and unsubstituted sydnones are almost equal (see Supplementary Information) which supports synchronous cyclization mechanism. At the same time, the N1-C11 bond in sydnones 10-12 (X = Cl, Br, I) is shorter than the C3-C10 bond. This supports our hypothesis that steric clashes can lower reaction rate in case of bulky substituent X.

Figure 7 .

 7 Figure 7. Transition state of cycloaddition reaction 12, between sydnone (X=I) and BCN. Only few atoms of both reactants are shown for clarity. Dashed line shows steric contacts H…I between halogen atom of sydnone and H atoms attached to one of carbons at the triple bond of alkyne, dotted lines indicate new single bonds forming in this reaction.

Table 1 .

 1 Rate constants of reaction shown on Figure1. The measurements were performed in water -DMSO mixture.

	No	R	X	Cycloalkyne

Rate constant (tolerance), M - 1 *sec -1

	1	p-MeO C 6 H 4	H	BCN	0.006 (0.001)	[14]
	2	p-Me C 6 H 4	H	BCN	0.032 (0.001)	[15]
	3	C 6 H 5	H	BCN	0.027 (0.002)	[14]
	4	p-CO 2 H C 6 H 4 H	BCN	0.059 (0.001)	[14]
	5	p-CF 3 C 6 H 4	H	BCN	0.199 (0.002)	[14]
	6	p-NO 2 C 6 H 4	H	BCN	0.289 (0.012)	[14]
	7	C 6 H 5	CH 3	BCN	0.018 (0.002)	[14]
	8	C 6 H 5	C 6 H 5	BCN	0.027 (0.001)	[14]
	9	C 6 H 5	CF 3	BCN	0.008 (0.001)	[14]
	10	p-Me C 6 H 4	Cl	BCN	0.872 (0.034)	[15]
	11	p-Me C 6 H 4	Br	BCN	0.592 (0.021)	[15]
	12	p-Me C 6 H 4	I	BCN	0.306 (0.008)	[15]
	13	p-CO 2 H C 6 H 4 Br	BCN	0.798 (0.065)	[14]
	14	p-CO 2 H C 6 H 4 Cl	BCN	1.593 (0.034)	[14]
	15	p-Me C 6 H 4	F	BCN	42	[15]
	16	p-Me C 6 H 4	F	TMTH	1500	[15]
	17	p-F C 6 H 4	F	TMTH	3500	[15]
	18	p-CF 3 C 6 H 4	F	TMTH	12000	[15]
						Reference

  where : 6 @ -Hirshfield charge on the atom A in molecule with M electrons, ;number of electrons in neutral molecule. Geometry of cation-radical and anionradical molecules corresponded to lowest energy structure of neutral molecule.

			010 . 2310 /2
			-. 2310 . /010
	The HOMO and LUMO energies were obtained within DFT scheme [20] for
	sydnone molecules.
	Fukui indexes are used to characterize atom ability to share/withdraw
	electronic charge. For calculation of Fukui nucleophilicity 5 6 7 , elecrophilicity 5 6 8
	and radical attack susceptibility 5 6 9 indexes of atom A single point calculations of
	molecules with added and removed electron were done:
		5 6 7	: 6 ;	: 6 ; 1
		5 6 8	: 6 ; 1	: 6 ;
	5 6 9	0.5 * : 6 ; 1	: 6 ; 1

Acknowledgement

The authors thank the Agence Nationale de la Recherche (ANR ClickReal, 2014-2019) for support.


Introduction to the special section:

When migrants produce the city: Everyday negotiations of urban space Danièle Bélanger, Nadine Cattan, Diego Coletto, Jean-Baptiste Frétigny This special section analyzes the dynamics and relations between migrants and city dwellers with the aim of moving beyond paradigms that focus on migrant victimization, exploitation, and precarity. At the crossroad of migration studies and urban studies, this collection seeks to explore how migrants and informal workers strategically negotiate their lives in the city and how they become an important part of the urban fabric. This section pays special attention to the informal aspects of how these social actors negotiate relationships and challenges in inhospitable environments. It brings together research on spatial and social relations between migrants and other actors in the public spaces of global metropolises, as well as secondary cities, and highlights the multiple ways in which the urban official order is challenged and subverted on an everyday basis. This collection focuses on migrants' and informal workers' agency and the contributions they make to the construction of urban life in Paris, Lima, Izmir (Turkey), Mexico City, and Tijuana.

Researchers have tended to overlook the theoretical, methodological, and empirical benefits of working at the intersection of migration studies and urban studies. The six articles in this section open up new paths by bridging these fields. All six manuscripts consider the rights to the city and to mobility for people in informal situations, such as being undocumented, having an ambiguous migratory status, and/or working in the informal sector. This special section seeks to address the following points. How do migrants and other city dwellers confront and redefine their identity in diverse public places? What kind of everyday strategies do they mobilize in order to enhance their access to the city and make advantageous arrangements on multiple scales? Which social networks improve the appropriation of public space and how is this space used in everyday activities? What conclusions can be drawn concerning the contested imaginaries of the city as these are negotiated, (re)produced, and challenged by migrants, actors in the informal economy, residents, tourists, urban planners, and governments, both local and national? A key hypothesis explored in this special section is that these capabilities, social relations, and manifestations are modes of political engagement with both symbolic and material features of the city.

The authors in this special section go beyond questioning how well accepted migrants are in cities to explain how sedentary and mobile regimes are linked in global metropolises and analyze the place attributed-or not attributed-to persons on the move.

The notion of temporariness and transit as modes of mobility are unsettled by the empirical findings of the articles. They show how populations on the move, while disrupting the regulation and imaginary of cities, participate in urban life, challenging modes of governance that are profoundly transformed by the intensification of all types of movement. The emotional content of this debate relates to considerations mainly of security, which challenge the right to the city by the legitimate and illegitimate orders of different actors. Each article shows how migrants, exiled persons, and refugees all claim a right to the city and highlights how these groups act on a city's materiality by questioning the sedentary urban moral order where transient individuals are seen as disruptive elements whose right to intervene in the production of the city must be denied.

Guillermo Candiz and Tanya Basok study the strategies developed by asylum seekers from Central America who are waiting in Tijuana, the border town with the United States. They use the research paradigm of mobile commons to highlight the dilemma faced by migrants in the trade-offs to be made in this border zone. They are caught between two choices: asking for access to the city and its resources by claiming citizenship rights, which give them a certain degree of protection (linked notably to obtaining migratory status which is recognized in Mexico), or pursuing their migration project and continuing their journey to the United States, despite further obstacles, via what the authors call "acts of mobility." This tension in transit territories, of which Tijuana is emblematic, shows not only how migratory policies but also migrants' urban experiences transform such cities into hybrid spaces, where waiting is mingled with frontierization and so powers mobility as an active and ambivalent component of the urban fabric.

In their article about migrant street vendors in Paris (France), Nadine Cattan and Jean-Baptiste Frétigny also question the right to the city through actions; however, they do so by surveying these migrants' capacities to manage interpersonal spaces with other users of urban spaces, especially tourists. By analyzing the emotional experiences stemming from interactions between vendors and tourists, they highlight three types of acceptability linked to bargaining, exoticism, and charity. Their article shows how, within this encounter of two types of globalization, the migrants contribute to reconnecting emblematic tourist spaces in Paris with the social diversity of a metropolis in all its political and ethical dimensions. Street vendors' creativity is demonstrated in their mastery of the distance and proximity they obtain through ephemeral access to the city.

The issue of access to urban resources by migrants is examined in very different circumstances by Laurent Faret, Andrea Paula González Cornejo, Jéssica Natalia Nájera Aguirre, and Itzel Abril Tinoco González. The migrants studied are en route to the United States and are often trapped in Mexico City, which is more and more at the confluence of a wide variety of migrant flows. This article focuses on the daily practices of newly arrived migrants from Central America and the Caribbean. It shows how the increasingly restrictive migration policy in Mexico City leads migrants to develop alternative, often informal, strategies in order to gain access to housing, work, health care, and education. Yet, even when they obtain elementary resources, they remain very vulnerable. The authors stress that such claims are all part of shared urban experiences, calling for a re-examination of the right to the city as a strategic and political instrument. Diego Coletto, Lema Jaber, and Linus Vanhellemont examine rights to the city in terms of the unprecedented lockdown situations brought about by the COVID-19 pandemic. By studying two markets in Lima (Peru), they reveal the reactivity and adaptability to this substantial change by street vendors, some of whom are migrants from rural provinces. These two food markets, serving poorer populations, continued to operate during lockdown periods. Sales activities allowed vendors to acquire a new visibility and enter into a renewed dialogue with the local authorities about the relocation of their activities in Lima's public spaces. The article, thus, goes beyond the victimization image of these informal actors to pursue a new understanding and de-homogenization of their situation. By highlighting the importance of difference, it is possible to better ascertain the issues involved in this urban fabric and the public policies that seek to address a great variety of informal situations in which migrants are often stakeholders.

The territorial scale of the deployment of informal urban fabric stands out in Annaelle Piva's article, which examines the spatiality in camp settlements of exiled persons in Paris. She highlights the power games at work in setting up the limits of non-rights to the city and its resources. By making camp-sites invisible, police action represents a twofold regulation of space through dispersion and "peripherization." Despite the life-cycle of camp settlements reflecting the transitory nature and vulnerability of refugees, the author, drawing on the theories of Rancière, reveals alternative strategies for appropriating urban space and the emancipation processes refugees engage in to claim an effective space in towns. The article shows how camps become transgressive places that challenge the order established by the authorities and how they reflect the political dimension of exiled persons' presence in Paris's urban spaces.

By relocating the debate on the governance of refugees from urban margins (camps) to city centres, Cenk Saraçoğlu and Danièle Bélanger study the co-presence of Syrian refugees and Turkish citizens in Izmir, a major Turkish coastal city on the Aegean. This study aligns with Cattan and Frétigny's work as it also observes that migrants assert their presence in these urban centres. It examines the daily life of Syrians in poorer neighbourhoods, characterized by a large variety of actors. An impressive body of ethnographical material makes it possible to decipher two forms of regulation on Syrians' presence in towns, established by host communities and resident Turks, namely socio-spatial marginalization and physical violence. Using a detailed, road-by-road scale of the city's districts, the article provides evidence of how control is established to discipline Syrian refugees' daily lives, by maintaining them in a powerfully unequal situation regarding access to urban resources. The way in which the issue of the legitimacy and acceptability of informal behaviour in the urban fabric is regulated leads the authors to conclude that there has been a depoliticization of the urban society formed by refugees and the host communities.

These six articles go beyond the binaries that dominate discussions about the lives of migrants and informal workers by shedding light on the hybridization processes at play in the city. Original insights into the everyday urban lives are offered, with an emphasis on both the (in)visibility of migrants and people who live in informal situations and the (in)formal relationships they construct with other actors in public spaces. These articles use a time-space perspective that attends to everyday (im)mobilities and itineraries. They document the ways in which different categories of informal workers and migrants access public spaces and, by (in)formally negotiating their right to the city, are both potential and effective co-producers of urban places and societies.