Paloma Vidal-Matutano
email: paloma.vidal@ulpgc.es

Auréade Henry

Yolanda Carrión Marco

Ethel Allué

P Vidal

Disentangling human from natural factors: Taphonomical value of microanatomical features on archaeological wood and charcoal assemblages

Archaeobotanical charcoal and wood analyses rely on the observation of different macro and microanatomical features affecting wood structure to variable extents. These features may result from a wide range of intrinsic and extrinsic factors alluding to different stages of the wood's taphonomical history: initial growth conditions, human selection, transformation/use and discard, post-depositional processes and archaeological sampling strategies. Papers in this volume address taphonomy in this broad sense, through recent methodological work mainly based on experimentation and case studies from a variety of chrono-cultural and geographical contexts. The authors present a number of tools available to wood and charcoal analysts and discuss their archaeological relevance to characterize anthropogenic and/or natural processes. The presented approaches are complementary and well reflect the extent to which wood and charcoal remains provide new insights into past human practices and social dynamics.

Introduction

Since the beginning of systematic charcoal analysis or anthracology, increasing attention has been given to different types of microanatomical features to advance the discipline's theoretical, methodological and interpretative framework [START_REF] Vernet | Les charbons de bois, les anciens écosystèmes et le rôle de l'homme[END_REF][START_REF] Thiébault | Charcoal analysis: methodological approaches, palaeoecological results and wood uses: proceedings of the Second International Meeting of Anthracology[END_REF][START_REF] Damblon | Proceedings of the Fourth International Meefing of Anthracology: Brussels[END_REF]; Fiorentino and Magri, 2008;Badal et al., 2012;[START_REF] Ludemann | Anthracology: local to global significance of charcoal science[END_REF][START_REF] Théry-Parisot | Anthracology and taphonomy, from wood gathering to charcoal analysis. A review of the taphonomic processes modifying charcoal assemblages, in archaeological contexts[END_REF][START_REF] Dotte-Sarout | Not just carbon: assessment and prospects for the application of anthracology in Oceania[END_REF][START_REF] Kabukcu | Wood Charcoal Analysis in Archaeology[END_REF]. In the early stages of anthracology, research was focused on sampling and quantification methods and potential representativeness biases such as fragmentation and mass reduction, which led to establishing the methodological basis of the discipline (Chabal, 1982[START_REF] Chabal | L'étude paléoécologique de sites protohistoriques à partir des charbons de bois: la question de l'unité de mesure[END_REF], 1992[START_REF] Chabal | Apports récents de l'anthracologie à la connaissance des paysages passés: performances et limites[END_REF], 1997;[START_REF] Chabal | L'Anthracologie[END_REF][START_REF] Krauss-Marguet | Analyse anthracologique du gisement Postglaciaire de la Poujade (Commune de Millau, Aveyron)[END_REF][START_REF] Badal García | Méthodes utilisées pour l'étude de sites préhistoriques[END_REF][START_REF] Badal | L'anthracologie préhistorique: à propos de certains problèmes méthodologiques[END_REF][START_REF] Pernaud | Paléoenvironnements végétaux et sociétés à l'Holocène dans le nord du Bassin Parisien: anthracoanalyses de sites archéologiques d'Île-de-France et de Picardie: méthodologie et paléoécologie[END_REF]. The establishment of an accurate method for charcoal analysis allowed to compare data between sites and to establish reliable paleoecological records [START_REF] Asouti | Reconstructing woodland vegetation and its exploitation by past societies, based on the analysis and interpretation of archaeological wood charcoal macroremains[END_REF]. Beyond these aspects, analyses of microanatomical features on charcoal developed in two main directions. First, the study of anatomical features to reconstruct firewood management, especially in Paleolithic and Mesolithic contexts for which a systemic interpreta tive framework was needed [START_REF] Théry-Parisot | Économie des combustibles au Paléolithique: expérimentation, taphonomie, anthracologie[END_REF]Henry and Théry-Parisot 2014a;[START_REF] Vidal-Matutano | Dead wood gathering among Neanderthal groups: charcoal evidence from Abric del Pastor and El Salt (Eastern Iberia)[END_REF]. Secondly, the study of anatomical features to identify aspects related to woodland management and paleoecological reconstructions [START_REF] Terral | Reconstruction of Holocene climate in southern France and eastern Spain using quantitative anatomy of olive wood and archaeological charcoal[END_REF][START_REF] Marguerie | Charcoal analysis and dendrology: data from archaeological sites in north-western France[END_REF][START_REF] Ludemann | Die Wälder am Schauinsland und ihre Nutzung durch Bergbau und Köhlerei[END_REF][START_REF] Dufraisse | Mesurer les diamètres du bois de feu en anthracologie. Outils dendrométriques et interprétation des données[END_REF][START_REF] Dufraisse | Anthraco-typology as a key approach to past firewood exploitation and woodland management reconstructions. Dendrological reference dataset modelling with dendro-anthracological tools[END_REF].

Macro-and microanatomical features visible within the ligneous cellular structure depend on different intrinsic and extrinsic factors alluding to various stages of the wood's taphonomical history: initial growth conditions, human selection, transformation/use and discard, post-depositional processes and archaeological sampling strategies. Initial growth conditions can be reflected by distinguishable anatomical variations of the wood structure resulting from the action of intrinsic or extrinsic natural factors or human activities such as deforestation, pruning or domestication/irrigation [START_REF] Schweingruber | Wood Structure and Environment[END_REF][START_REF] Dufraisse | Charcoal analysis: new analytical tools and methods for archaeology[END_REF]Terral, 2000;Limier et al., 2018). The identification of the features related to these processes sheds light on woodland management practices and environmental effects on woody plants. Human selection is a major topic in anthracology where the study of anatomical features allows going beyond the taxonomic identification to analyze charcoal assemblages, as aspects related to socio-economical behavior can be approached through studies aiming at reconstruct ing minimal diameters and state of the wood (green, seasoned, decayed) used. Indeed, the latter are additional indicators of firewood acquisition modalities (woodcutting vs. gathering) and have different implications in terms of heat requirements and fire use (Moskal-del [START_REF] Hoyo | Preservation of fungi in archaeological charcoal[END_REF][START_REF] Théry-Parisot | Seasoned or green? Radial cracks analysis as a method for identifying the use of green wood as fuel in archaeological charcoal[END_REF]Henry and Théry-Parisot, 2014b;Mallol and Henry, 2017;[START_REF] Vidal-Matutano | Firewood and hearths: Middle Palaeolithic woody taxa distribution from El Salt, stratigraphic unit Xb (Eastern Iberia)[END_REF]. Transformation/use and discard involve different processes subjected to variable archaeological visibility according to the studied material; macro-and microtraces on worked wood allow reconstructing manufacture processes and uses, whereas in the case of charcoal assemblages the combustion constitutes an additional taphonomy filter. Burning induces wood fragmentation and mass reduction, alters the cell walls and can produce specific microanatomical alterations such as cracks and vitrification, in relation to the initial fuel characteristics [START_REF] Chabal | L'Anthracologie[END_REF][START_REF] Théry-Parisot | Anthracology and taphonomy, from wood gathering to charcoal analysis. A review of the taphonomic processes modifying charcoal assemblages, in archaeological contexts[END_REF][START_REF] Théry-Parisot | Seasoned or green? Radial cracks analysis as a method for identifying the use of green wood as fuel in archaeological charcoal[END_REF][START_REF] Mcparland | Is vitrification in charcoal a result of high temperature burning of wood?[END_REF]. Post-depositional processes and preservation of charcoal or wood in archaeological or natural contexts are diverse. In the course of their depositional or post-depositional history, charcoal particles can be fragmented or hyper fragmented. Anthracological remains are affected by processes indicating the action of mechanical, biological or climatic factors [START_REF] Chrzazvez | Impact of post-depositional processes on charcoal fragmentation and archaeobotanical implications: experimental approach combining charcoal analysis and biomechanics[END_REF]. In the case of wood, its preservation generally implies specific depositional and post-depositional conditions. Finally, archaeological sampling strategies and recovery methods can have a significant impact on assemblage fragmentation and, more generally, representativeness, which stresses the importance of the traceability of excavation and laboratory methods [START_REF] Vidal-Matutano | Firewood and hearths: Middle Palaeolithic woody taxa distribution from El Salt, stratigraphic unit Xb (Eastern Iberia)[END_REF][START_REF] Arranz-Otaegui | Evaluating the impact of water flotation and the state of the wood in archaeological wood charcoal remains: implications for the reconstruction of past vegetation and identification of firewood gathering strategies at Tell Qarassa North (south Syria)[END_REF].

To understand processes and classify different features several methods have been applied relying overall on experimentation and ethnoarchaeology (Ntinou, 2002;Théry-Parisot et al., 2009;[START_REF] Henry | Ethnoarchaeology of fuel use in northern forests: towards a better characterization of prehistoric fire-related activities[END_REF]Théry Parisot et al., 2016;[START_REF] Gelabert | The ethnoarchaeology of firewood management in the Fang villages of Equatorial Guinea, central Africa: Implications for the interpretation of wood fuel remains from archaeological sites[END_REF][START_REF] Dussol | Experimental anthracology: evaluating the role of combustion processes in the representivity of archaeological charcoal records in tropical forests, a case study from the Maya Lowlands[END_REF][START_REF] Caruso-Fermé | Tracing Driftwood in Archaeological Contexts: Experimental Data and Anthracological Studies at the Orejas De Burro 1 Site (Patagonia, Argentina)[END_REF][START_REF] Mcparland | Is vitrification in charcoal a result of high temperature burning of wood?[END_REF][START_REF] Dufraisse | Anthraco-typology as a key approach to past firewood exploitation and woodland management reconstructions. Dendrological reference dataset modelling with dendro-anthracological tools[END_REF]. Accurate observations of the materials and efforts in understanding the processes involved in the formation of microanatomical features have provided interpretation keys of anthracological materials that are currently at different stages of development.

For the Palaeolithic and, more generally, temporary hunter-gatherer occupations, specific questionings related to taphonomy issues and wood-related human behavior in terms of fire management systems are widely being developed. These questionings concern mainly fuel selection and acquisition modalities from an environmental, territorial and functional perspective [START_REF] Théry-Parisot | Économie des combustibles au Paléolithique: expérimentation, taphonomie, anthracologie[END_REF]2002;[START_REF] Théry-Parisot | Anthracology and taphonomy, from wood gathering to charcoal analysis. A review of the taphonomic processes modifying charcoal assemblages, in archaeological contexts[END_REF], 2018;[START_REF] Allué | Preliminary issues regarding the taphonomic study of archaeological charcoal upon the record from Abric Romani (Capellades, España)[END_REF][START_REF] Allué | Charcoal taphonomy: the study of the cell structure and surface deformations of Pinus sylvestris type for the understanding of formation processes of archaeological charcoal assemblages[END_REF]Henry and Théry-Parisot, 2014a;[START_REF] Vidal-Matutano | Firewood and hearths: Middle Palaeolithic woody taxa distribution from El Salt, stratigraphic unit Xb (Eastern Iberia)[END_REF][START_REF] Vidal-Matutano | Neanderthal firewood management: evidence from stratigraphic unit IV of Abric Del Pastor (eastern Iberia)[END_REF][START_REF] Vidal-Matutano | The anthropogenic use of firewood during the European Middle Pleistocene: charcoal evidence from levels XIII and XI of Bolomor Cave, Eastern Iberia (230-160 ka)[END_REF][START_REF] Alcolea-Gracia | Donde hubo fuego: estudio de la gestión humana de la madera como recurso en el valle del Ebro entre el Tardiglaciar y el Holoceno Medio[END_REF]Allué & Mas, this volume;Martínez-Varea et al., this volume). These specific chrono-cultural contexts have also led researchers to investigate more in-depth the impact of depositional and post-depositional processes in order to contribute to reconstructing on-site taphonomic conditions (Théry-Parisot, 2001; [START_REF] Chrzazvez | Impact of post-depositional processes on charcoal fragmentation and archaeobotanical implications: experimental approach combining charcoal analysis and biomechanics[END_REF][START_REF] Marquer | Étude des macro-, méso-et micro-charbons du site épigravettien de Mezhyrich (Ukraine): données taphonomiques et anthracologiques[END_REF][START_REF] Vidal-Matutano | Dead wood gathering among Neanderthal groups: charcoal evidence from Abric del Pastor and El Salt (Eastern Iberia)[END_REF][START_REF] Vidal-Matutano | Firewood and hearths: Middle Palaeolithic woody taxa distribution from El Salt, stratigraphic unit Xb (Eastern Iberia)[END_REF]. With increasing anthropogenic pressure from the Neolithic onwards and within more recent, generally better preserved contexts, the development of dendroanthracological tools and multi-proxy approaches have allowed significant advances in the knowledge of past woodland management and agricultural systems [START_REF] Ludemann | Die Wälder am Schauinsland und ihre Nutzung durch Bergbau und Köhlerei[END_REF][START_REF] Carrión-Marco | Dendrología y arqueología: las huellas del clima y de la explotación humana de la madera[END_REF][START_REF] Delhon | Environment and landscape management during the Middle Neolithic in Southern France: evidence for agro-sylvo-pastoral systems in the Middle Rhone Valley[END_REF][START_REF] Dufraisse | Charcoal analysis: new analytical tools and methods for archaeology[END_REF][START_REF] Terral | Bio-archaeological evidence of olive tree (Olea europaea L.) irrigation during the Middle Ages in Southern France and North Eastern Spain[END_REF][START_REF] Delhon | Shepherds and plants in the Alps: multi-proxy archaeobotanical analysis of Neolithic dung from "La Grande Rivoire" (Isère, France)[END_REF][START_REF] Paradis-Grenouillet | Deciduous oak/chestnut: Differential shrinkage of wood during charcoalification? Preliminary experimental results and implications for wood diameter study in anthracology[END_REF][START_REF] Marguerie | Dendroanthracologie[END_REF][START_REF] Dufraisse | Mesurer les diamètres du bois de feu en anthracologie. Outils dendrométriques et interprétation des données[END_REF][START_REF] Dufraisse | Anthraco-typology as a key approach to past firewood exploitation and woodland management reconstructions. Dendrological reference dataset modelling with dendro-anthracological tools[END_REF][START_REF] Out | Branch age and diameter: useful criteria for recognising woodland management in the present and past?[END_REF]. Finally, exceptional contexts such as waterlogged sites where non-charred ligneous materials are preserved allow the adoption of a wider perspec tive regarding wood uses, from techno-functional and cultural approaches of ligneous materials and artefacts to the reconstruction of woodland management and economies [START_REF] Mertens | Hölzerne Fischfanggeräte und ihre Bedeutung für die Ökonomie des Mesolithikums[END_REF][START_REF] Chabal | Le mobilier organique des puits antiques et autres contextes humides de Lattara[END_REF][START_REF] Carrión | Revealing Iberian woodcraft: conserved wooden artefacts from south-east Spain[END_REF][START_REF] Klooß | They were fishing in the sea and coppicing the forest. Terminal Mesolithic and Early Neolithic wooden artefacts of coastal settlements on the southwestern Baltic Sea[END_REF][START_REF] López-Bultó | Wood procurement at the early Neolithic site of La Draga (Banyoles, Barcelona)[END_REF]Martín-Seijó and Carrión, 2012;[START_REF] Zhilin | Early art in the Urals: New research on the wooden sculpture from Shigir[END_REF].

Case studies presented in this special issue

The case studies presented in this volume cover a wide geographic range from Eurasia to South America. The archaeological contexts are also diverse, with a chronological span from the Middle Palaeolithic to historical periods. There is a broad application of new approaches to traditional anthracology that definitely contribute to a better understanding of anatomical features and, more generally, of the anthracological record. These approaches are intended to solve particular issues regarding specific human activities related to wood uses and to understand natural and anthropogenic processes affecting anthracological assemblages.

As previously mentioned, the Paleolithic has traditionally been the most studied period from a taphonomical perspective, probably because many processes related to charcoal aging and decay make evident the need for a deep taphonomic analysis. In this volume, two papers focus on this period from different perspectives. Martínez-Varea et al. present two case studies, Abrigo de la Quebrada and Cova de les Cendres (Spain) to analyze the taphonomical characteristics of charcoal remains affected by post-depositional processes. The authors analyze preservation conditions involving different type of contexts and human occupation patterns. Allué and Mas present an overview of microanatomical features observed on Pinus sylvestris type charcoal from Paleolithic sequences in Iberia. As the values provided by the microanatomical features are difficult to interpret in terms of human behavior, they are used as a comparative proxy between archaeological sequences. However, the advance in the description of these features and an accurate classification will allow in the future a better understanding of the meaning of these results with broader implications for wood uses during the Paleolithic.

Charcoal quantification and floristic representativeness of the anthracological spectra (floristic diversity and proportions between taxa) is an important issue in anthracology. The main variables identified as affecting the reliability of charcoal assemblages are the recovery methods used during the excavation, the number of analyzed fragments, the study of different contexts and the origin of the samples (dispersed vs concentrated charcoal assemblages) [START_REF] Chabal | L'étude paléoécologique de sites protohistoriques à partir des charbons de bois: la question de l'unité de mesure[END_REF]1992;1997;[START_REF] Chabal | L'Anthracologie[END_REF][START_REF] Badal | L'anthracologie préhistorique: à propos de certains problèmes méthodologiques[END_REF]. The general validity of these premises, mainly built around protohistorical to historical case studies, is reexamined by Henry et al. for a Mesolithic cave occupation. The objective was to understand the respective contributions of different size classes (2-4 mm and > 4 mm) and sampling contexts (dispersed vs. concentrated charcoal) on the floristic results of Escabasses Cave (France). The fact that the floristic results vary more according to size class than to sampling context goes against previous observations and invites taking fragment size more systematically into account. Furthermore, this specific situation may have paleothnographic implications regarding the way in which small mobile Mesolithic groups managed fireplaces in a confined space.

Combustion processes entailing wood carbonization have major implications in different aspects affecting archaeological charcoal assemblages. Some work has already been done using experimental archaeology to understand specific processes such as mass reduction, wood anatomy changes and fragmentation (Slocum et al. 1978;[START_REF] Chrzazvez | Impact of post-depositional processes on charcoal fragmentation and archaeobotanical implications: experimental approach combining charcoal analysis and biomechanics[END_REF][START_REF] Gonçalves | The effect of carbonization on wood structure of Dalbergia violacea, Stryphnodendron polyphyllum, Tapirira guianensis, Vochysia tucanorum, and Pouteria torta from the Brazilian cerrado[END_REF]Théry-Parisot et al., 2010b). Among the alterations resulting from combustion, vitrification is one of the most recurrent, even though it rarely affects the entire assemblage. This glassy and refractive aspect of the wood anatomical cell struc ture has no accurate explanation yet, even though most of the studies agree that this alteration results frombut cannot be only explained by -the combustion process. Different parameters have been tested and analyzed by several researchers, such as type of combustion process, temperature or initial state of the wood [START_REF] Théry-Parisot | Économie des combustibles au Paléolithique: expérimentation, taphonomie, anthracologie[END_REF][START_REF] Mcparland | Is vitrification in charcoal a result of high temperature burning of wood?[END_REF][START_REF] Henry | Paléoenvironnements et gestion du bois de feu au Mésolithique dans le sud-ouest de la France : anthracologie, ethnoarchéologie et expérimentation[END_REF]. [START_REF] Courty | Fuel origin and firing product preservation in archaeological occupation contexts[END_REF] proposed that this alteration results from wood flash-pyrolysis induced by aeroplasma. In this volume, Courty et al. present research focused on vitrification affecting different kinds of charcoal remains. The authors present new evidence based on modern burning assemblages and their identification and a detailed process inducing vitrified charcoal. One of the strengths of the study is to provide insights on fuel qualities and pyrotechnology, showing that there are different degrees, and probably different processes, leading to the production of vitrified charcoal components. Vidal Matutano et al. present a multi-approach study based on the spatial analysis of charcoal and burned faunal remains to understand the causes of vitrification in an anthracological assemblage from the pre-Hispanic site of Chasogo (Tenerife, Spain). In this study, abundant vitrified wood charcoal fragments seem to be spatially related to the combustion structures and the burnt bones suggesting the relationship of this phenomenon with specific combustion conditions, allowing the discussion of some working hypotheses for future analyses.

Wood management has become an important issue based overall in dendroanthracology, with diameter estimates used for drawing paleoeconomical but also palaeoenvironmental hypotheses regarding woodland structure [START_REF] Carrión-Marco | Dendrología y arqueología: las huellas del clima y de la explotación humana de la madera[END_REF][START_REF] Dufraisse | Mesurer les diamètres du bois de feu en anthracologie. Outils dendrométriques et interprétation des données[END_REF]. Dufraisse et al. present new methodological developments for the analysis of wood diameter estimation using the Pith Location Tool (PLT). This method, based on trigonometry, was applied to seven different taxa, constituting an important reference dataset. The strengths of this research are to provide measurement protocols and guidelines for an accurate dendroanthracological analysis, as well as an R-function for wood and charcoal fragments allowing a wide use of this tool among the research community.

Caruso-Fermé and Théry-Parisot apply for the first time an evaluation of the minimal wood caliber to a Patagonian hunter-gatherer site, Cerro Casa de Piedra 7 (Argentina). They propose a comparison between diameters measured on charcoal vs. wood fragments of Nothofagus, the most abundant taxon in this archaeological context. The results show completely different patterns between firewood and wood diameters used at the site throughout time, which has interesting implications regarding wood acquisition strategies of these Patagonian societies.

López-Bultó presents the famous Neolithic waterlogged wood assemblage of La Draga (Spain), from which he analyses woodworking debris and formless pieces of wood. Indeed, in sites with exceptional preservational conditions, this last category is often left aside in xylological studies. This study confirms the necessity to take them into account, as they allow the discussion of key aspects of wood economies such as raw material selection, transport and manufacturing processes (e.g., Caruso Fermé 2012; [START_REF] Taylor | The Wooden Artefacts[END_REF].

Finally, post-depositional processes can be studied from a wide perspective including pedoanthracology as an approach to understand soil formation processes in natural contexts [START_REF] Scott | Charcoal: Taphonomy and significance in geology, botany and archaeology[END_REF]. Ponomarenko et al. present a multidisciplinary study based on sedimentology, palynology and pedoanthracology. The study is focused on the discrimination of natural versus anthropogenic causes of erosion phases during historical times in the Middle Volga region. According to the data, different land uses and cultivation techniques are clearly recognizable in the sedimentary profiles of the investigated area. The authors conclude that human practices related to livestock and deforestation are the main cause of erosion.

Concluding thoughts

Anthracology is a relatively young research field that is still expanding thanks to new analytical methods mainly based on experimental work and actualistic reference datasets in general. Ongoing approaches range from the macroscopic to the molecular levels, with a whole new world of possibilities opening up for a better understanding of both natural and anthropogenic processes. The articles in this SI only represent partially what is going on out there, but in our opinion, they well reflect the fact that charcoal and wood analyses tend to diversify around multiple topics. As shown in this volume, these topics go from paleoenvironmental reconstructions and land-use patterns to very specific questions regarding firewood and wood management, while methodological developments and questionings are still at the heart of the discipline. Much is still left to do to be able to characterize, understand and interpret the many signatures identifiable within the charcoal record, and their relevance for the study of different taxa, archaeological levels, sites and periods. Taphonomic issues are often treated laterally in papers focused on other topics, which often makes it difficult to follow the progress made. In this sense, this volume aims to be a meeting point for the definition, evaluation of archaeological relevance and general discussion of microanatomical signatures observed in archaeological charcoal and wood. [START_REF] Badal | Análisis antracológico de los restos del fuego doméstico del abrigo de los Baños (Ariño, Teruel)[END_REF], [START_REF] Badal | La presencia de hongos e insectos xilófagos en el carbón arqueológico: propuestas de interpretación[END_REF]

Uncited references

  , Caruso-Ferme and Thery-Parisot (this volume), Courty et al. (this volume), Dufraisse et al. (this volume), Henry et al. (this volume), Lopez-Bulto (this volume), Ponomarenko et al. (this volume), Vidal-Matutano et al. (this volume).

Acknowledgements

We would like to thank the European Association of Archaeologists for agreeing to host the session 110 "Disentangling human from natural factors: taphonomical value of microanatomical features on archaeological wood and charcoal assemblages" at the EAA Annual Meeting 2018 in Barcelona. Many thanks to the presenters of the 12 oral communications and the poster, as well as to the colleagues who attended the session. Thanks to A. Aşınmaz for his volunteer work as a technician in our session room.

Special thanks to the research projects and institutions that have financially supported the organization of the session:

GDR 3544 "Sciences du bois" and CEPAM-CNRS UMR 7264 laboratory (A. Henry and P. Vidal-Matutano); research project MINECO Ramón y Cajal program RYC-2014-16027 (Y. Carrión Marco); research projects MICINN-FEDER PGC2018-093925-B-C32, MINECO (HAR2016-76760-C3-1-P) and the Generalitat de Catalunya research project 2017 SGR 836 (E. Allué).

Many thanks to JAS: Reports for agreeing to host this Special Issue about charcoal taphonomy, and especially to Chris O. Hunt and Vinaya Sundar.

We also would like to express our deepest gratitude to our colleagues and peers, the 17 reviewers who put a lot of effort in providing useful and pertinent comments, which contributed to improve all of the papers composing this volume. It is important to acknowledge their work, as they agree to perform these evaluation tasks without compensation of any kind, often within an already overloaded work schedule. Thank you for your support, as this Special Issue would not have been possible without your hard work: Cláudia Costa (ICArEHB-FCHS, University of Algarve, Portugal), Claire Delhon (Université Côte d'Azur, CNRS, CEPAM UMR 7264, France), Thomas Feiss (SILVA laboratory, UMR 1434, University of Lorraine, France), Tryfon Giagkoulis (Aristo


Daniele Loco

Riccardo Spezia
email: riccardo.spezia@sorbonne-universite.fr

François Cartier

Isabelle Chataigner
email: isabelle.chataigner@univ-rouen.fr

Jean-Philip Piquemal
email: jean-philip.piquemal@sorbonne-universite.fr

Solvation Effects Drive the Selectivity in Diels-Alder Reaction Under Hyperbaric Conditions

come   L'archive ouverte pluridisciplinaire

The favourable effect of high pressure (HP) regimes on the kinetics of certain organic reactions is known for a while, 1 notably for hindered substrates. HP indeed activates reactions unfeasible under conventional conditions, often increasing chemo-, regioand stereo-selectivities. [2][START_REF] Van Eldik | High Pressure Chemistry: Synthetic, Mechanistic, and Supercritical Applications[END_REF][START_REF] Chataigner | Maddaluno in Activation methods: Sonochemistry and High Pressure[END_REF] It also induces a significant reduction in energy costs by circumventing the continuous heating of the medium. As such, it can be considered as a "physical catalysis".

HP effect has been markedly studied on the Diels-Alder cycloaddition, [START_REF] Matsumoto | [END_REF]6 a remarkable chemical transformation which still provides key solutions for the synthesis of (poly)cyclic molecules. The inertia of the components may limit its applicability, for steric or electronic reasons. The process is then usually activated by heating, using catalysts or changing solvents. Temperature increase often favors the reverse reaction, while catalysts may induce substrates and/or products degradation. Working under pressure is often an efficient alternative.

Cycloaddition reactions are characterized by very negative activation volumes (∆V ‡ = -20 to -40 cm [START_REF] Van Eldik | High Pressure Chemistry: Synthetic, Mechanistic, and Supercritical Applications[END_REF] •mol -1 ) 7 and their kinetics are usually accelerated under high pressure. One of the first and most emblematic examples is the use of furan as a diene in the synthesis of cantharidin by Dauben (see Figure 1, top). 8,9 Since then, HP activation has shown its efficiency in performing (4+2) cycloadditions involving sterically hindered [10][11][12] or electronically deactivated [13][14][15][16] substrates for instance. This soft and eco-compatible activation mode has thus given access to different complex and functionalized structures, in the context of total syntheses for example. 17,18 The positive influence of pressure on the diastereoselectivity of Diels-Alder cycloaddition reactions has been observed in many cases (see Figure 1, bottom). Generally, the more compact endo approach is favored by a HPs. 2,19,20 In some cases endo kinetic adducts are obtained under hyperbaric cycloaddition, while exo thermodynamic ones are generated by heating. 21 The activation observed in organic reactions in liquid phase is explained by the Evans-Polanyi equation, which connects a negative activation volume with the acceleration of the reaction with pressure increase. Based on such a theory, Ladanyi and Hynes studied pressure effects on rate constant of simple systems, finding an important dependence of the activation volume at very HP. 22 Later, Dumas et al. have studied an aza-Michael addition via an atomistic model where the effect of pressure was included modifying the van der Waals radii. 23 More recently, Hoffmann, Cammi and coworkers have studied the effect of pressure from a quantum chemistry approach coupled with an implicit solvation model on a set of organic reactions. 24 The activation volumes was thus evaluated from electronic densities.

To get more insight into this phenonenon, we undertook a computational modeling study of pressure effects on a prototypical Diels-Alder reaction, between cyclopentadiene and acrylonitrile (see Figure 2). Since the HP effect can modify the solute, the solvent and their interactions, we have adopted a physically grounded description of the solvent by means of an advanced classical polarizable force field (FF), 26 coupled with a quantum chemistry method to treat the solute reacting species, which are then explicitly affected by the presence of the classical but atomistic solvent. This is what we define a polarizable variational QM/MM method. [27][28][29] To sample the large number of solute-solvent configurations needed to extract meaningful energies and properties from such atomistic models we resorted to a sequential approach: first performing full FF based molecular dynamics to have equilibrated solvent configurations at the needed temperature and high/low pressures, then including the quantum treatment of the solute in the structure sampling and performing shorter QM/MM MDs, starting from different points of the full FF-based trajectories.

The solute molecules are described within the the density functional theory (DFT) approximation, with the M06-2X functional 30 and a 3-21G atomic basis set. A detailed discussion for this choice is reported in the Supporting Informations (SI, Section S1). Note that it is not uncommon that a small basis set can provide results better than a larger one: 31 convergence in basis set sometimes occur only for a very large one, which would be then inadequate to perform QM/MM MDs.

During the polarizable MDs the solvent has been modeled with the AMOEBA FF, which performs correctly at different pressures without any ad hoc parametrization. 32 Solute was solvated by 603 CH 2 Cl 2 molecules using an AMOEBA potential recently developed for organic molecules. 33 The solvent was first equilibrated at room temperature and at two different pressures (low, 1 atm and high, 10 katm) fixing the solutes in their QM optimized geometry and running about 2 ns polarizable MD simulations using the Bussi-Parrinello algorithm to keep an NPT ensemble. 34 This procedure was performed for each chemical species relevant to describe the reaction pathway, and thus placed as solute inside the CH 2 Cl 2 box of solvent. These structures represent: (i) the reactants, (ii) transition states (TS), endo and exo approaches and (iii) cycloadducts (CA), endo and exo diastereomers. Intermediates between reactants and TS were also considered (vide infra) and named here (iv) pre-TS. These structures were taken from the intrinsic reaction coordinate (IRC) which connects reactants with TS (see SI for more details). The solvent equilibration was performed for the two different pressure regimes.

After an equilibration time of less than 10 ps the systems reached the equilibrium density at the two different pressures, namely 1.419(3) and 1.776(4) g/mL. These values are independent on the chemical nature of the solutes. In Figure 3 the density time series at the beginning of the classical polarizable trajectories are reported, showing the differences between low and high pressure simulations, for which also different simulation box sizes were obtained. Finally, the solute-solvent interaction was refined via a variational formulation of hybrid polarizable QM/MM, accounting for the mutual polarization between the classical and the quantum subsystems. The forces needed for QM/MM dynamics were consistently derived from the energies computed in this way. [27][28][29]35 This approach was used for the first time here to investigate pressure effects on a chemical reaction.

Once the solvation was equilibrated via polar-izable MD, QM/MM simulations were performed, either running few picoseconds QM/MM MD or single point QM/MM energy calculations (see Section S3 of SI), using structures sampled from the full-FF based polarizable MDs. The Tinker-HP package 36 coupled with a development version of Gaussian09 was used for all the calculations. [START_REF] Frisch | Gaussian Development Version[END_REF] The reaction energy profiles were thus obtained extracting average energies from the simulations at different pressures, setting at zero the energy of the reactants at both 1 atm and 10 katm. Such relative energies are reported in Figure 4 where we also draw the structures of important geometries along the reaction pathway.

In the 1 atm reaction pathway, the TS was found to be 11 and 18 kcal/mol higher than reactants for endo and exo approaches, respectively (with a statistical uncertainty of 1 kcal/mol for these and further energy values). Products were exothermic of 38 and 39 kcal/mol for endo and exo pathways. Note that for such a reaction a very small difference between the two stereomers is expected. [START_REF] Larranaga | [END_REF] Finally, at 1 atm the pre-TS has an energy lying in between the reactants and the TS, as in the IRC calculations from which the structures were extracted.

At 10 katm the situation is very different instead: the TS was found to be energetically more stable than the reactants. Interestingly, the pre-TS structures are now more stable than the corresponding TS. However, while for the exo reaction, a barrier of about 9 kcal/mol was calculated between the pre-TS and the TS, for the endo reaction it was found to be almost zero (1 kcal/mol which lies into the uncertainty). This means that the rate constant in the case of the endo approach is almost under diffusion control and much faster than the exo one. Finally, products were similar in energy with a slight preference for the endo diastereomer.

Results show that the reaction is indeed kinetically favoured by a pressure increase. For the exo approach, an energy barrier of 18 kcal/mol at ambient pressure translates into a 9 kcal/mol one at HP; for endo approach, the difference is even bigger: the 11 kcal/mol barrier observed at ambient pressure almost vanishes at HP and the kinetics is then diffusion controlled. Finally, the resulting endo cycloadduct is favoured (due to kinetic control) with respect to exo, as experimentally expected.

In our microscopic description of the system, increasing the pressure triggers a compact packing of the molecules. As a consequence, the attractive interactions with the solvent molecules increase (see S2 of SI for the absolute energetics). In this respect, the use of polarizable FFs 26 is crucial: the polarization variation with the pressure ensures a large part of this stabilization (see Table S2 of SI). Only thanks to the fluctuations of their induced dipole, the classical solvent molecules are allowed to respond to the different pressure regimes and readjust during the compaction of the molecular packing.

On one hand, the extended exo-stereomers are more exposed to the solvent molecules than the compact endo ones. This effect is even greater in the pre-TS where the substrates are further apart: the stabilizing interactions induced by the solvent under HP show their maximum effect for the exo pre-TS, which is more stable than the endo one by 2 kcal/mol. On the other hand the TS endo is generally considered to be more favoured than the TS exo, due to the attractive secondary orbital interactions between the p orbitals that are not directly involved in the formation of the new σ bonds, according to the famous (even if still debated) proposition by Woodward and Hoffman 39 . As a net result of these two antagonistic phenomena, the endo TS is almost barrier less, while the exo one features an energy activation of 9 kcal/mol. Thus, the solvation effect explains the greater preference for the endo cycloadduct formation under HP.

We connect our results to the commonly accepted picture offered by the concept of activation volume. Following the approach used by Cammi and co-workers, 40 we compute molar volumes from QM/MM electron densities at the two different pressures, using them to compute activation and reaction volumes. As shown in Table 1, for each species the activation volume does not depend neither on the pressure, nor on the stereochemistry. This means that we are in the linear regime of pressure dependence, confirming experimental data. [START_REF] Van Eldik | High Pressure Chemistry: Synthetic, Mechanistic, and Supercritical Applications[END_REF] Furthermore, the reaction volumes are similar, showing the high similarity between TS and products as expected for such a late transition state Diels-Alder reaction.

The pre-TS were found to be crucial in differenciating the kinetics between the exo and endo approaches at HP. This is in line with the different activation volumes calculated when the pre-TSs were considered as reactants (see Table 1). The exo pre-TS volume is indeed similar to the reactants one, while the endo pre-TS volume resembles more the corresponding TS at HP (the absolute volume of each structure is reported in Fig. S4 of SI).

In conclusion, we propose here that the solvation plays a fundamental role on the pressure effect on this Diels-Alder prototypical cycloaddition. Under HP, the efficient solvation of the TSs lowers the energy barriers and explain the acceleration of the process. In addition to this phenomenon, the exo pre-TS arrangement, more exposed to solvent molecules, is favorably solvated under HP and more stabilized than the endo one, making the passage to the corresponding TS more challenging. In the meantime, the compact endo TS remains favored thanks to its attractive secondary orbital interactions. These two parallel phenomena explain the enhancement of the diastereoselectivity and the greater endo preference for Diels-Alder cycloaddition under HP. Note that this interpretation is somehow in line with the work by Houk and Bickelhaupt, showing that a distortion/interaction model is especially fruitful for explaining cycloaddition reactivity. 41 The QM/MM method developed in our group shows how pressure effects in the condensed phase can be modeled for an archetypal Diels-Alder reaction. The original approach proposed will certainly find other applications, not only for transformations under HP, but also to further investigate the crucial role of solvents/additives in organic reactions.

Figure 1 :

 1 Figure 1: Examples of successful selectivity enhancement in organic chemistry reactions.

25

 25

Figure 2 :

 2 Figure 2: Model Diels-Alder reaction used in the present work. TS: transition state, CA: cycloadduct

Figure 3 :

 3 Figure 3: Top) Solution density's time-series along the initial part of the full-FF equilibration dynamics run for pre-TSs, TSs, cycloadducts and reactants in exo configuration, in CH2Cl2 solution. In red are reported the densities for the 10 katm MDs, while those for 1 atm MDs are in blue. Bottom) Shrinking of the MD simulation box due to the pressure increase. The cubic boxes reported are taken from two trajectories of the same solvated system at 1 and 10 katm, with a box side length of 39 and 36 Å respectively

Figure 4 :

 4 Figure 4: Top) Structures of the different molecular species involved in the reaction. Bottom) Energy diagram as from QM/MM simulations

Table 1 :

 1 Activation volume ∆V ‡ and reaction volume ∆V react for the reaction under investigation computed from the SCF QM/MM electron density with Monte-Carlo integration by Gaussian 0940 .

	The values,

Acknowledgements

This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 810367), project EMC2. Computations have been performed at GENCI on the Occigen machine (CINES, Montpellier, France) on grant no A0070707671.

Keywords


domain: sdv phys

text:
Anastasia A Minervina

Ekaterina A Komech

Aleksei Titov

Bensouda Meriem

Elisa Koraichi

Ilgar Rosati

Andre Mamedov

Grigory Franke

Dmitriy Efimov

Thierry Chudakov

Mora

Meriem Bensouda Koraichi

Elisa Rosati

Ilgar Z Mamedov

Andre Franke

Grigory A Efimov

Dmitriy M Chudakov

Thierry Mora

Aleksandra M Walczak

Yuri B Lebedev

Mikhail V Pogorelyy

Longitudinal high-throughput TCR repertoire profiling reveals the dynamics of T cell memory formation after mild COVID-19 infection

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

INTRODUCTION

COVID-19 is a global pandemic caused by the novel SARS-CoV-2 betacoronavirus [START_REF] Vabret | Immunology of COVID-19: current state of the science[END_REF]. T cells are crucial for clearing respiratory viral infections and providing longterm immune memory [START_REF] Schmidt | The CD8 T Cell Response to Respiratory Virus Infections[END_REF][START_REF] Swain | Expanding roles for CD4+ T cells in immunity to viruses[END_REF]. Two major subsets of T cells participate in the immune response to viral infection in different ways: activated CD8+ T cells directly kill infected cells, while subpopulations of CD4+ T cells produce signaling molecules that regulate myeloid cell behaviour, drive and support CD8 response and the formation of long-term CD8 memory, and participate in the selection and affinity maturation of antigen specific Bcells, ultimately leading to the generation of neutralizing antibodies. In SARS-1 survivors, antigen-specific memory T cells were detected up to 11 years after the initial infection, when viral-specific antibodies were undetectable [START_REF] Ng | Memory T cell responses targeting the SARS coronavirus persist up to 11 years post-infection[END_REF][START_REF] Oh | Engineering T Cells Specific for a Dominant Severe Acute Respiratory Syndrome Coronavirus CD8 T Cell Epitope[END_REF]. The T cell response was shown to be critical for protection in SARS-1-infected mice [START_REF] Zhao | T Cell Responses Are Required for Protection from Clinical Disease and for Virus Clearance in Severe Acute Respiratory Syndrome Coronavirus-Infected Mice[END_REF]. Patients with X-linked agammaglobulinemia, a genetic disorder associated with lack of B cells, have been reported to recover from symptomatic COVID-19 [START_REF] Quinti | A possible role for B cells in COVID-19? Lesson from patients with agammaglobulinemia[END_REF][START_REF] Soresina | Two X-linked agammaglobulinemia patients develop pneumonia as COVID-19 manifestation but recover[END_REF], suggesting that in some cases T cells are sufficient for viral clearance. Theravajan et al. showed that activated CD8+HLA-DR+CD38+ T cells in a mild case of COVID-19 significantly expand following symptom onset, reaching their peak frequency of 12% of CD8+ T cells on day 9 after symptom onset, and contract thereafter [START_REF] Thevarajan | Breadth of concomitant immune responses prior to patient recovery: a case report of non-severe COVID-19[END_REF]. Given the average time of 5 days from infection to the onset of symptoms [START_REF] Bi | Epidemiology and transmission of COVID-19 in 391 cases and 1286 of their close contacts in Shenzhen, China: a retrospective cohort study[END_REF], the dynamics and magnitude of T cell response to SARS-CoV-2 is similar to that observed after immunization with live vaccines [START_REF] Miller | Human effector and memory CD8+ T cell responses to smallpox and yellow fever vaccines[END_REF]. The exact immunodominant CD8+ and CD4+ SARS-CoV-2 epitopes are yet unknown. However, SARS-CoV-2-specific T cells were detected in COVID-19 survivors by activation following stimulation with SARS-CoV-2 proteins [START_REF] Ni | Detection of SARS-CoV-2-Specific Humoral and Cellular Immunity in COVID-19 Convalescent Individuals[END_REF], or by viral protein-derived peptide pools [START_REF] Weiskopf | Phenotype of SARS-CoV-2specific T-cells in COVID-19 patients with acute respiratory distress syndrome[END_REF][START_REF] Braun | Presence of SARS-CoV-2 reactive T cells in COVID-19 patients and healthy donors[END_REF]. Some of the T cells activated by peptide stimulation were shown to have a memory phenotype [START_REF] Weiskopf | Phenotype of SARS-CoV-2specific T-cells in COVID-19 patients with acute respiratory distress syndrome[END_REF], and some potentially cross-reactive CD4+ T cells were found in healthy donors [START_REF] Braun | Presence of SARS-CoV-2 reactive T cells in COVID-19 patients and healthy donors[END_REF][START_REF] Grifoni | Targets of T cell responses to SARS-CoV-2 coronavirus in humans with COVID-19 disease and unexposed individuals[END_REF]. T cells recognise short pathogen-derived peptides presented on the cell surface of the Major Histocompatibility Complex (MHC) using hypervariable T cell receptors (TCR). TCR repertoire sequencing allows for the quantitative tracking of T cell clones in time, as they go through the expansion and contraction phases of the response. It was previously shown that quantitative longitudinal TCR sequencing is able to identify antigen-specific expanding and contracting T cells in response to yellow fever vaccination with high sensitivity and specificity [START_REF] Minervina | Primary and secondary antiviral response captured by the dynamics and phenotype of individual T cell clones[END_REF][START_REF] Pogorelyy | Precise tracking of vaccineresponding T cell clones reveals convergent and personalized response in identical twins[END_REF][START_REF] Dewitt | Dynamics of the cytotoxic T cell response to a model of acute viral infection[END_REF]. Not only clonal expansion but also significant contraction from the peak of the response are distinctive traits of T cell clones specific to the virus [START_REF] Pogorelyy | Precise tracking of vaccineresponding T cell clones reveals convergent and personalized response in identical twins[END_REF].

In this study we use longitudinal TCRalpha and TCRbeta repertoire sequencing to quantitatively track T cell clones that significantly expand and contract after a mild COVID-19, and determine their phenotype. We reveal the dynamics and the phenotype of the memory cells formed after infection, identify pre-existing T cell memory clones participating in the response, and describe public TCR sequence motifs of SARS-nCoV-2-reactive clones, suggesting a response to immunodominant epitopes.

RESULTS

In the middle of March (day 0) donor W female and donor M (male, both healthy young adults), returned to their home country from the one of the centers of the COVID-19 outbreak in Europe at the time. Upon arrival, according to local regulations, they were put into strict self-quarantine for 14 days. On day 3 of self-isolation both developed low grade fever, fatigue and myalgia, which lasted 4 days and was followed by a temporary loss of smell for donor M. On days 15, 30, 37 and 45 we collected peripheral blood samples from both donors. The presence of SARS-CoV-2 specific antibodies in the plasma was measured at all timepoints using SARS-CoV-2 S-RBD domain specific ELISA (Fig. S1). From each blood sample we isolated PBMCs (peripheral blood mononuclear cells, in two biological replicates), CD4+, and CD8+ T cells. Additionally, on days 30 and 45 we isolated four T cell memory subpopulations (Fig. S2): Effector Memory (EM: CCR7-CD45RA-), Effector Memory with CD45RA reexpression (EMRA: CCR7-CD45RA+), Central Memory (CM: CCR7+CD45RA-), and Stem Cell-like Memory (SCM: CCR7+CD45RA+CD95+). From all samples we isolated RNA and performed TCRalpha and TCRbeta repertoire sequencing as previously described [START_REF] Pogorelyy | Persisting fetal clonotypes influence the structure and overlap of adult human T cell receptor repertoires[END_REF]. For both donors, TCRalpha and TCRbeta repertoires were obtained for other projects one and two years prior to infection. Additionally, TCR repertoires of multiple samples for donor M -including sorted memory subpopulations -are available from a published longitudinal TCR sequencing study after yellow fever vaccination (donor M1 samples in [START_REF] Minervina | Primary and secondary antiviral response captured by the dynamics and phenotype of individual T cell clones[END_REF]).

From previously described activated T cell dynamics for SARS-CoV-2 [START_REF] Thevarajan | Breadth of concomitant immune responses prior to patient recovery: a case report of non-severe COVID-19[END_REF], and immunization with live vaccines [START_REF] Miller | Human effector and memory CD8+ T cell responses to smallpox and yellow fever vaccines[END_REF], the peak of the T cell expansion is expected around day 15 post-infection, and responding T cells significantly contract by day 45. However, Weiskopf et al. [START_REF] Weiskopf | Phenotype of SARS-CoV-2specific T-cells in COVID-19 patients with acute respiratory distress syndrome[END_REF] reports an increase of SARS-CoV-2-reactive T cells at later timepoints, peaking in some donors after 30 days following symptom onset. To identify groups of T cell clones with similar dynamics in an unbiased way, we used Principal Component Analysis (PCA) in the space of T cell clonal trajectories (Fig. 1b andc). In both donors, and in both TCRalpha and TCRbeta repertoires, we identified three clusters of clones with distinct dynamics. The first cluster (Fig. 1bc, purple) corresponded to abundant TCR clonotypes which had constant concentrations across timepoints, the second cluster (Fig. 1bc, green) showed contraction dynamics from day 15 to day 45, and the third cluster (Fig. 1bc, yellow), showed an unexpected clonal expansion from day 15 with a peak on day 37 followed by contraction. The clustering and dynamics are similar in both donors and are reproduced in TCRbeta (Fig. 1bc) and TCRalpha (Fig. S3) repertoires. We next used edgeR, a software for differential gene expression analysis [START_REF] Robinson | edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[END_REF], to specifically detect changes in clonotype concentration between pairs of timepoints in a statistically reliable way. EdgeR uses biological replicate samples collected at each timepoint to train a negative-binomial noise model. We identified 374 TCRalpha and 373 TCRbeta clonotypes in donor W, and 797 TCRalpha and 672 TCRbeta in donor M significantly contracted from day 15 to day 45 (largely overlapping with cluster 2 of clonal trajectories). 438 TCRalpha and 533 TCRbeta for donor W, and 531 TCRalpha and 599 TCRbeta clonotypes for donor M were significantly expanded from day 15 to 37 (corresponding to cluster 3 of clonal trajectories). Reproducing the analysis using NoiseET, a Bayesian differential expansion model [START_REF] Touzel | Inferring the immune response from repertoire sequencing[END_REF], yielded similar results (see Fig. S4), suggesting that our results are robust to the choice of statistical model and parameter inference procedure. Note that, to identify putatively SARS-CoV-2 reactive clones, we only used post-infection timepoints, so that our analysis can be reproduced in other patients and studies where pre-infection timepoints are unavailable. However, tracking the identified responding clones back to pre-infection timepoints reveals strong clonal expansions from pre-to post-infection (Fig. 1de, Fig. S3cd). For brevity, we further refer to clonotypes significantly contracted from day 15 to 45 as contracting clones and clonotypes significantly expanding from day 15 to 37 as expanding clones. For each contracting and expanding clone we determined their CD4/CD8 phenotype using separately sequenced repertoires of CD4+ and CD8+ subpopulations (see Methods). Both CD4+ and CD8+ subsets participated actively in the response (Fig. 1de). Interestingly, clonotypes expanding after day 15 were significantly biased towards the CD4+ phenotype, while contracting clones had balanced CD4/CD8 phenotype fractions in both donors (Fisher exact test, p < 0.0001 for both donors).

On days 30 and 45 we identified both contracting (Fig. 2a,b) and expanding (SI Fig. 5a-c) T cell clones in the memory subpopulations of peripheral blood. Both CD4+ and CD8+ responding clones were found in the CM and EM subsets, however CD4+ were more biased towards CM, and CD8+ clones much more represented in the EMRA subset. A small number of both CD4+ and CD8+ responding clonotypes were also identified in the SCM subpopulation, which was previously shown to be a longlived T cell memory subset [START_REF] Marraco | Long-lasting stem celllike memory CD8+ T cells with a naïve-like profile upon yellow fever vaccination[END_REF]. Intriguingly, a number of responding CD4+ clones, but very few CD8+ clones, were also represented in the repertoires of both donors 1 and 2 years before the infection. Pre-existing clones were expanded after infection, and contracted afterwards for both donors (SI Fig. 6). For donor M, for whom we had previously sequenced memory subpopulations before the infection [START_REF] Minervina | Primary and secondary antiviral response captured by the dynamics and phenotype of individual T cell clones[END_REF], we were able to identify pre-existing SARS-CoV-2-reactive CD4+ clones in the CM subpopulation 1 year before the infection. Interestingly, on day 30 after infection the majority of pre-infection CM clones were detected in the EM subpopulation, suggesting recent T cell activation and a switch of the phenotype from memory to effector. These clones might represent memory T cells cross-reactive for other infections, e.g. other human coronaviruses. The search for responding clones TCRbeta amino acid sequences in VDJdb [START_REF] Bagaev | VDJdb in 2019: database extension, new analysis infrastructure and a T-cell receptor motif compendium[END_REF] -a database of TCRs with known specificity -resulted in essentially no overlap: a total of 3 matches all corresponded to the EBV (Epstein-Barr virus) epitope presented by the HLA-A*03 MHC allele, which is absent in both donors (SI Table 1). At the time of the analysis there were no SARS-1, SARS- CoV-2 or any other coronavirus-specific TCR sequences in VDJdb, so the absence of matches suggests that contracting and expanding clones are unlikely to be specific for immunodominant epitopes of common pathogens covered in VDJdb.

It was previously shown that TCRs recognising the same antigens frequently have highly similar TCR sequences [START_REF] Dash | Quantifiable predictive features define epitope-specific T cell receptor repertoires[END_REF][START_REF] Glanville | Identifying specificity groups in the T cell receptor repertoire[END_REF]. To identify motifs in TCR amino acid sequences, we plotted similarity networks for significantly contracted (Fig. 2d-g) and expanded (Fig. S5d-g) clonotypes. In both donors we found clusters of highly similar clones in both CD4+ and CD8+ subsets for expanding and contracting clonotypes. Clusters were largely donor-specific, as expected, since our donors have dissimilar HLA alleles (SI Table 1) and thus each is likely to present a non-overlapping set of T cell antigens. The largest cluster, described by the motif TRAV35-CAGXNYGGSQGNLIF-TRAJ42, was identified in donor M's CD4+ contracting alpha chains. Clones from this cluster constituted 15.3% of all of donor M's CD4+ responding cells on day 15, suggesting a response to an immunodominant CD4+ epitope in the SARS-CoV-2 proteome. The high similarity of the TCR sequences of responding clones in this cluster allowed us to independently identify motifs from donor M's CD4 alpha contracting clones using the ALICE algorithm [START_REF] Pogorelyy | Detecting T cell receptors involved in immune responses from single repertoire snapshots[END_REF] (Fig. S8). While the time dependent methods (Fig. 1) identify abundant clones, the ALICE approach is complementary to both edgeR and NoisET as it identifies clusters of T cells with similar sequences independently of their individual abundances. Finding highly similar responding TCR sequences in both donors, despite their different HLA types, suggests that some of the identified responding clones could be public, and are likely to be found in other SARS-CoV-2 infected individuals. Using a mechanistic model of TCR recombination [START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF][START_REF] Sethna | Population variability in the generation and thymic selection of T-cell repertoires[END_REF] we predicted the probability of occurrence of the reactive clones (see Methods). For TCRbeta we also looked for the identified contracting and expanding clones in a large public dataset of TCRbeta repertoires from [START_REF] Emerson | Immunosequencing identifies signatures of cytomegalovirus exposure history and HLAmediated effects on the T cell repertoire[END_REF]. The TCRbeta occurrence probability correlated well with the frequency of contracting and expanding clones in this public dataset (Fig. S7), suggesting that at least a fraction of the SARS-CoV-2 responding clones are public and are likely to be shared between patients.

DISCUSSION

Using longitudinal repertoire sequencing, we identified a group of CD4+ and CD8+ T clones that contract after recovery from a coronavirus infection. Our response timelines agree with T cell dynamics reported by Theravajan et al. [START_REF] Thevarajan | Breadth of concomitant immune responses prior to patient recovery: a case report of non-severe COVID-19[END_REF] for mild COVID-19, as well as with dynamics of T cell response to live vaccines [START_REF] Miller | Human effector and memory CD8+ T cell responses to smallpox and yellow fever vaccines[END_REF]. Surprisingly, in both donors we also identified a group of predominantly CD4+ clonotypes which expanded from day 15 to day 37 after the infection. One possible explanation for this second wave of expansion is the priming of CD4+ T cells by antigen-specific B-cells, but there might be other mechanisms such as the migration of SARS-CoV-2 specific T cells from lymphoid organs. It is also possible that later expanding T cells are triggered by another infection, simultaneously and asymptomatically occurring in both donors around day 30. The major limitation of our study is that we do not know the specificity of contracted and expanded clones observed after infection. However, accumulation of TCR sequences from SARS-CoV-2-specific T cells will allow us to connect the TCR motifs we described to particular epitopes of the virus and estimate their relative contribution to the immune response. Databases of SARS-CoV-2-associated TCR sequences combined with sequencing of patients' TCR repertoires will allow for the tracking of the adaptive immune response in the early stages of the disease, and possibly identify features correlating with clinical outcomes. We showed that a large fraction of putatively SARS-CoV-2 reactive T cell clones are later found in memory subpopulations, and a subset of CD4+ clones were identified in pre-infection central memory subsets. The presence of SARS-CoV-2 crossreactive CD4 T cells in healthy individuals was recently demonstrated [START_REF] Braun | Presence of SARS-CoV-2 reactive T cells in COVID-19 patients and healthy donors[END_REF][START_REF] Grifoni | Targets of T cell responses to SARS-CoV-2 coronavirus in humans with COVID-19 disease and unexposed individuals[END_REF]. Our data further suggests that cross-reactive CD4+ T cells can participate in the response in vivo. It is interesting to ask if the presence of cross-reactive T cells before infection is linked to the mildness of the disease. Larger studies with cohorts of severe and mild cases with pre-infection timepoints are needed to address this question.

METHODS

Donors and blood samples

Peripheral blood samples from two young healthy adult volunteers, donor W (female) and donor M (male) were collected with written informed consent in a certified diagnostics laboratory. Both donors gave written informed consent to participate in the study under the declaration of Helsinki. HLA alleles of both donors (SI Table 1) were determined by an in-house cDNA high-throughput sequencing method.

SARS-CoV-2 S-RBD domain specific ELISA

An ELISA assay kit developed by the National Research Centre for Hematology was used for detection of anti-S-RBD IgG according to the manufacturer's protocol. The relative IgG level was calculated by dividing the OD (optical density) values by the mean OD value of the cut-off positive control serum supplied with the Kit. OD values of 2 samples (d37 and d45) for donor M exceeded the limit of linearity for the Kit. In order to properly compare the relative IgG levels between d30, d37 and d45, these samples were diluted 1:400 instead of 1:100, the ratios d37:d30 and d45:d30 were calculated and used to calculate the relative IgG level of d37 and d45.

Isolation of PBMCs and T cell subpopulations

PBMCs were isolated with the Ficoll-Paque density gradient centrifugation protocol. CD4+ and CD8+ T cells were isolated from PBMCs with Dynabeads CD4+ and CD8+ positive selection kits respectively. For isolation of EM, EMRA, CM and SCM memory subpopulations we stained PBMCs with the following antibody mix: anti-CD3-FITC (UCHT1, eBioscience), anti-CD45RA-eFluor450 (HI100, eBioscience), anti-CCR7-APC (3D12, eBioscience), anti-CD95-PE (DX2, eBioscience). Cell sorting was performed on FACS Aria III, all four isolated subpopulations were lysed with Trizol reagent immediately after sorting.

TCR library preparation and sequencing

TCRalpha and TCRbeta cDNA libraries preparation was performed as previously described in [START_REF] Pogorelyy | Persisting fetal clonotypes influence the structure and overlap of adult human T cell receptor repertoires[END_REF]. RNA was isolated from each sample using Trizol reagent according to the manufacturer's instructions. A universal primer binding site, sample barcode and unique molecular identifier (UMI) sequences were introduced using the 5'RACE technology with TCRalpha and TCRbeta constant segment specific primers for cDNA synthesis. cDNA libraries were amplified in two PCR steps, with introduction of the second sample barcode and Illumina TruSeq adapter sequences at the second PCR step. Libraries were sequenced using the Illumina NovaSeq platform (2x150bp read length).

TCR repertoire data analysis Raw data preprocessing. Raw sequencing data was demultiplexed and UMI guided consensuses were built using migec v.1.2.7 [START_REF] Shugay | Towards error-free profiling of immune repertoires[END_REF]. Resulting UMI consensuses were aligned to V and J genomic templates of the TRA and TRB locus and assembled into clonotypes with mixcr v.2.1.11 [START_REF] Bolotin | MiXCR: software for comprehensive adaptive immunity profiling[END_REF]. See SI Table 2 for the number of cells, UMIs and unique clonotypes for each sample.

Identification of clonotypes with active dynamics. Principal component analysis (PCA) of clonal trajectories was performed as described before [START_REF] Minervina | Primary and secondary antiviral response captured by the dynamics and phenotype of individual T cell clones[END_REF]. First we selected clones which were present among the top 1000 abundant in any of post-infection PBMC repertoires. Next, for each clone we calculated its frequency at each post-infection timepoint and divided this frequency by the maximum frequency of this clone for normalization. Then we performed PCA on the resulting normalized clonal trajectory matrix and identified three clusters of trajectories with hierarchical clustering with complete linkage, using Euclidean distances between trajectories. We identify statistically significant contractions and expansions with edgeR as previously described [START_REF] Pogorelyy | Precise tracking of vaccineresponding T cell clones reveals convergent and personalized response in identical twins[END_REF], using FDR adjusted p < 0.01 and log 2 fold change threshold of 1. NoisET implements the Bayesian detection method described in [START_REF] Touzel | Inferring the immune response from repertoire sequencing[END_REF]. Briefly, a two-step noise model accounting for cell sampling and expression noise is inferred from replicates, and a second model of expansion is learned from the two timepoints to be compared. The procedure outputs the posterior probability of expansion or contraction, and the median estimated log 2 fold change, whose thresholds are set to 0.05 and 1 respectively.

TCR generative probability estimation. Models of V(D)J-recombination for the TCRalpha and TCRbeta loci were inferred for each individual using IGoR [START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF]. Selection models were then inferred and evaluated using SONIA [START_REF] Sethna | Population variability in the generation and thymic selection of T-cell repertoires[END_REF] to estimate TCR occurrence probabilities.

FIG. 1 .

 1 FIG. 1.Longitudinal tracking of T cell clones after mild COVID-19. a, Study design. Peripheral blood of two donors was sampled longitudinally on days 15, 30, 37, 45 after arrival in Russia. At each timepoint, we evaluated SARS-CoV-2-specific antibodies using ELISA (Fig.S1) and isolated PBMCs in two biological replicates. Additionally, CD4+ and CD8+ cells were isolated from a separate portion of blood, and EM, CM, EMRA, SCM memory subpopulations were FACS sorted on days 30 and 45. For each sample we sequenced TCRalpha and TCRbeta repertoires. For both donors pre-infection PBMC repertoires were sampled in 2018 and 2019 for other projects. b,c, PCA of clonal temporal trajectories identifies three groups of clones with distinctive behaviours. Left: first two principal components of the 1000 most abundant TCRbeta clonotype normalized frequencies in PBMC at post-infection timepoints. Right: mean ± 2.96 SE of the normalized trajectory of TCR clones in each cluster. Contracting (d) and expanding (e) clones include both CD4+ and CD8+ T cells, and are less abundant in pre-infection repertoires. T cell clones significantly contracted from day 15 to day 45 (d) and significantly expanded from day 15 to day 37 (e) were identified using edgeR in both donors. The fraction of contracting (d) and expanding (e) TCRbeta clonotypes in the total repertoire (corresponding to the fraction of responding cells of all T cells) is plotted in log-scale for all reactive clones (left), reactive clones with the CD4 (middle) and CD8 (right) phenotypes. Similar dynamics were observed in TCRalpha repertoires (Fig.S3), and for significantly expanded/contracted clones identified with the NoisET Bayesian differential expansion statistical model (Fig.S4).

FIG. 2 .

 2 FIG.2. Memory phenotypes and TCR sequence motifs of responding clonotypes contracting after day 15. a, A large fraction of contracting clonotypes is identified in T cell memory subsets after infection. Fraction of contracting CD4+ and CD8+ TCRbeta clonotypes present in 1-year and 2-year pre-infection PBMC and in at least one of memory subpopulation sampled on day 30 and day 37 post infection. b, Responding clones are found in different memory subsets. For both W (left) and M (right) donors, CD4+ clonotypes were more frequently found in the Central Memory (CM) subset than CD8+ T cells, and both CD4+ and CD8+ contracting clonotypes were present in the Effector Memory (EM) compartment. c, For donor M (right), CD4+ (but few CD8+) contracting clonotypes are also identified in memory subsets 1 year before the infection, with a bias towards the CM subpopulation. d, Analysis of TCR amino acid sequences of contracting clones reveal distinctive motifs. Each vertex in the similarity network corresponds to a contracting clonotype. An edge indicates 2 or less amino acid mismatches in the CDR3 region (and identical V and J segments). Networks are plotted separately for CD4alpha (d), CD4beta (e), CD8alpha (f ) and CD8beta (g) contracting clonotypes. Clonotypes without neighbours are not shown. Sequence logos corresponding to the largest clusters are shown under the corresponding network plots.

FIG

  FIG. S1. Both donors developed anti-SARS-CoV-2 IgG response by day 30 post infection. The relative level of SARS-CoV-2 S-RBD domain specific IgG (y-axis) is plotted against time. Solid black line shows the threshold for positive testing.

  FIG. S7. TCR occcurence probability predicts public responding clonotypes. a Probability of TCR generation and selection (Ppost) of TCRbeta clones with active dynamics (both contracted and expanded from both donors) is strongly correlated with the occurrence of these clones in the large public database of TCRbeta repertoires from Emerson et al. [29]. (b)Distribution for TCR occurency probability of SARS-CoV-2 responding TCRbeta clonotypes. Dashed line shows TCR occurrence probability distribution for all clones in total repertoire. Some identified clones are public and have large recombination probabilities, and thus are expected to be shared with additional SARS-CoV-2 patients.

ACKNOWLEDGMENTS

The study was supported by 00351. E.R and A.F. are supported by the DFG Excellence Cluster Precision Medicine in Chronic Inflammation (Exc2167) and the DFG grant n. 4096610003, M.B.K., T.M., and A.M.W. are supported by European Research Council COG 724208. I.Z.M. is supported by RFBR 19-54-12011 and 18-29-09132. D.M.C. is supported by the grant from Ministry of Science and Higher Education of the Russian Federation 075-15-2019-1789.

Data availability

Raw sequencing data are deposited to the Short Read Archive (SRA) accession: PRJNA633317. Processed TCRalpha and TCRbeta repertoire datasets, resulting repertoires of SARS-CoV-2-reactive clones, and raw data preprocessing instructions can be accessed from: https: //github.com/pogorely/Minervina_COVID.


text:
Giulio Isacchini

Carlos Olivares

Armita Nourmohammad

Aleksandra M Walczak

Thierry Mora

SOS: Online probability estimation and generation of T and B cell receptors

 and SONIA [3]-have been publicly released to the community that allow for the inference of generative and selection models from high-throughput sequencing data. However using these tools requires some scripting or command-line skills and familiarity with complex datasets. As a result the application of the above models has not been available to a broad audience. In this application note we fill this gap by presenting Simple OLGA & SONIA (SOS), a web-based interface where users with no coding skills can compute the generation and post-selection probabilities of their sequences, as well as generate batches of synthetic sequences. The application also functions on mobile phones. Availability and implementation:

INTRODUCTION

The adaptive immune systems recognises pathogens through the generation of a highly diverse repertoire of T and B cell receptors (TCR and BCR) which have the potential to recognise even unknown pathogen and initiate an immune response. In order to produce this diversity it exploits a highly stochastic process named V(D)J recombination. In addition, to block possible auto-reactive receptors, a selection process is mounted in the thymus for T cells, and a similar process of central tolerance is implemented for B cells. Probabilistic models of TCR and BCR have been proposed [START_REF] Murugan | Statistical inference of the generation probability of Tcell receptors from sequence repertoires[END_REF][START_REF] Ralph | Consistency of VDJ Rearrangement and Substitution Parameters Enables Accurate B Cell Receptor Sequence Annotation[END_REF][START_REF] Elhanati | Quantifying selection in immune receptor repertoires[END_REF] based on immune repertoire sequencing data [START_REF] Georgiou | The promise and challenge of high-throughput sequencing of the antibody repertoire[END_REF][START_REF] Heather | Highthroughput sequencing of the T-cell receptor repertoire: pitfalls and opportunities[END_REF][START_REF] Minervina | TCR and BCR repertoire profiling in adaptive immunity[END_REF][START_REF] Bradley | Using T Cell Receptor Repertoires to Understand the Principles of Adaptive Immune Recognition[END_REF]. Software has been developped to infer the probability of generation of any B-or T-cell receptor (IGoR [START_REF] Marcou | Highthroughput immune repertoire analysis with IGoR[END_REF]), and to evaluate this probability for both nucleotide and amino-acid sequences (OLGA [START_REF] Sethna | OLGA: fast computation of generation probabilities of B-and T-cell receptor amino acid sequences and motifs[END_REF]). Another tool (SONIA [START_REF] Sethna | Population variability in the generation and thymic selection of T-cell repertoires[END_REF]) was released to infer the selective pressures acting on the receptors and used to predict the probability of naive sequences in the periphery [START_REF] Isacchini | On generative models of T-cell receptor sequences[END_REF]. In order to make these tools available to a broader audience, we provide a new web tool which allows for the analysis of single TCR and BCR sequences.

FEATURES

As explained in the introductory "About" tab, the web tool evaluates the generation and post-selection probability of single naive T and B cell receptors in different species based on the specific sequence the user inputs manually. The engine is based on two pieces of python software, OLGA and SONIA, and shipped with pre-trained models of recombination and selection for the following loci: human alpha and beta chains or TCR (TRA and TRB), human heavy and light chain of unmutated BCR (IGH, IGK, and IGL), and mouse TRB.

After choosing the species and receptor chain in the "Evaluate" tab, the user inputs a Complementary Determining Region 3 (CDR3), either as a nucleotide or an amino acid sequence, and optionally V and J germline genes from dropdown lists. The server outputs the generation probability (P gen , conditioned on sequence productivity), and the post selection probability (P post ), as shown in Fig. 1 (left). When V and J are not specified, the program sums over all possibilites for these segments to calculate the total probability of the CDR3.

To help interpret the result and assess how the sequence of interest compares to others, P gen , P post , and the selection factor Q = P post /P gen are plotted as green vertical lines on histograms of random sequences taken pre-(blue line) and post-(orange line) selection (Fig. 1, right). That feature only works when V and J and specified. The tool also provides an estimation of the probability to observe the sequence in a generic repertoire. The user inputs the size N of the sequenced repertoire (unique productive nucleotide sequences), and the tool outputs the probability of observing the sequence within a repertoire of that size, given by 1 -(1 -P post ) N .

Using the "Generate" tab, the user can synthetize a specified number of receptor sequences from P gen or P post , after choosing the species and chain type from dropdown lists. The file with the generated sequences, composed of the CDR3 sequence (nucleotide and aminoacid translation), V and J segments, is available for download as a CSV file. The user may fix the seed of the random number generator for reproducibility. arXiv:2003.13149v1 [q-bio.GN] 29 Mar 2020 FIG. 1: SOS web interface. The user inputs a CDR3 sequence (amino acid or nucleotides) and V and J segments. The program outputs the generation probability Pgen, the probability in the periphery Ppost, and evaluates a p-value corresponding to the probability of finding that sequence by chance in a repertoire of size N (input by user). An additional tab allows for the generation of synthetic repertoires.

DISCUSSION

The interface can be used by investigators to evaluate how surprised one should be to find a given sequence in one or multiple repertoires. It could help distinguish receptors with a specific function from chance detections. The tool can also be used to evaluate the potential of certain receptors (in particular antibodies, albeit in their unmutated version) for vaccination or therapeutic purposes. The web interface is also available on mobile phones without the plotting options.

Acknowledgments. This work was partially supported by the European Research Council Proof of Concept Grant n. 824735

SOS is freely available to use at sites.google.com/view/ statbiophysens/sos with source code at github.com/statbiophys/sos.