domain: chim phys sdv

text:
Phuong Anh Dang

Carla Palomino-Durand

Mohamed Elsafi Mabrouk

Pierre Marquaille

Clément Odier

Sophie Norvez

Emmanuel Pauthe

Laurent Corté

Rational formulation design of injectable thermosensitive chitosan-based hydrogels for cell encapsulation and delivery

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

ABSTRACT (150 words)

This work reports a rational design of injectable thermosensitive chitosan systems for cell encapsulation and delivery. Using mixtures of two phosphate salts, beta-glycerophosphate and ammonium hydrogen phosphate, we demonstrate that the pH and the osmolarity can be adjusted separately by varying the molar ratios between the salts and the D-glucosamine monomers. We found the existence of a critical temperature above which gelation time decays following a power-law. This gelation kinetics can be finely tuned through the pH and salt-glucosamine ratios.

Formulations having physiological pH and osmolarity were produced for chitosan concentrations ranging from 0.4 to 0.9 wt%. They remain liquid for more than 2 h at 20°C and form a macroporous gel within 2 min at 37°C. In vitro encapsulation of pre-osteoblastic cells and gingival fibroblasts showed homogeneous cell distribution and good cell viability up to 24 h.

Such an approach provides a valuable platform to design thermosensitive cell-laden systems.

INTRODUCTION

Cell-laden injectable hydrogels have raised a great interest during the last decades as valuable delivery vehicles or matrices for cell therapies [START_REF] Li | Injectable and biodegradable hydrogels: Gelation, biodegradation and biomedical applications[END_REF]. Administrated in a minimally invasive manner, these gels fill-in the defects in hardly accessible and fragile tissues, and might provide a temporary scaffold to injected cells to ensure their viability and stimulate their activity. For effective clinical applications, such systems should stay liquid before and during injection, then gel rapidly in physiological conditions for temperature, pH and osmolarity.

In the liquid state, pH and osmolarity should be finely controlled for the survival of the encapsulated cells. After gelation, the gel scaffold must provide appropriate mechanical and chemical properties as well as a macroporosity enabling cell invasion and growth. Different stimuli-responsive systems for in situ gelling have been investigated, triggered by temperature [START_REF] Yap | Thermo-reversible injectable hydrogel composing of pluronic F127 and carboxymethyl hexanoyl chitosan for cell-encapsulation[END_REF], pH (Z. [START_REF] Li | PH-Sensitive and Thermosensitive Hydrogels as Stem-Cell Carriers for Cardiac Therapy[END_REF], or light [START_REF] Neves | Injectable Hyaluronic Acid Hydrogels Enriched with Platelet Lysate as a Cryostable Off-the-Shelf System for Cell-Based Therapies[END_REF]. Among those, thermogelling hydrogels based on biocompatible chitosan have been extensively studied [START_REF] Mekhail | Injectable Chitosan-Based Scaffolds in Regenerative Medicine and their Clinical Translatability[END_REF].

Chitosan (CS) is a linear copolymer composed of D-glucosamine and N-acetyl-D-glucosamine monomers (Figure 1). CS is obtained by alkaline deacetylation of parent chitin, a primary component of cell walls in fungi and crustacean shells. As one of the few representatives of biodegradable and biocompatible polymers approved for clinical use [START_REF] Berger | Structure and interactions in covalently and ionically crosslinked chitosan hydrogels for biomedical applications[END_REF][START_REF] Riva | Chitosan and chitosan derivatives in drug delivery and tissue engineering[END_REF][START_REF] Saravanan | A review on injectable chitosan/beta glycerophosphate hydrogels for bone tissue regeneration[END_REF] CS has long been a popular biomaterial for developing cell delivery matrices (Di [START_REF] Martino | Chitosan: A versatile biopolymer for orthopaedic tissueengineering[END_REF][START_REF] Zhou | Glycerophosphate-based chitosan thermosensitive hydrogels and their biomedical applications[END_REF]. Amino groups of chitosan present an apparent pKaapp in the vicinity of 6.2. At pH < pKaapp, CS behaves as a cationic polyelectrolyte bearing protonated amino groups, soluble in aqueous media. Near neutral pHs where CS is mostly deprotonated, the strong tendency of the carbohydrate polymer to self-associate through H-bonding and hydrophobic interactions dominates, giving rise to the formation of a hydrated gel-like precipitate [START_REF] Rinaudo | Chitin and chitosan: Properties and applications[END_REF]. This poor solubility of CS at physiological pHs has long been a major limitation for biomedical uses, until the pioneering work by Chenite et al., who demonstrated that in presence of beta-glycerophosphate (βGP, Figure 1), cationic CS solutions turned into thermally sensitive solutions, remaining liquid at room temperature even at neutral pHs, and forming hydrogels around 37°C [START_REF] Chenite | Novel injectable neutral solutions of chitosan form biodegradable gels in situ[END_REF][START_REF] Chenite | Rheological characterisation of thermogelling chitosan/glycerol-phosphate solutions[END_REF]. The mechanisms responsible for the thermoresponsivity of CS-βGP mixtures are still debated.

Cho et al. associated the gelling mechanism with enhanced hydrophobic interactions with increasing temperatures [START_REF] Cho | Physical gelation of chitosan in the presence of βglycerophosphate: The effect of temperature[END_REF][START_REF] Rinaudo | Chitin and chitosan: Properties and applications[END_REF]. Filion et al. underlined that the pKaapp of CS decreases by heating, producing a release of protons that may be accepted by the negatively charged βGP, a weak base having a pKa,2 of 6.65 at 25°C. The CS neutralization induced by the temperature-dependent proton transfer would result in a phase separation of the chitosan solution [START_REF] Filion | Chitosan-glycerol-phosphate (GP) gels release freely diffusible GP and possess titratable fixed charge[END_REF][START_REF] Lavertu | Heat-induced transfer of protons from chitosan to glycerol phosphate produces chitosan precipitation and gelation[END_REF]. Thermodynamic insights into the mechanism of the phase transition of CS induced by βGP have been recently related by Grinberg et al. (2020). From calorimetric analyses, they concluded that the βGP binding induces the formation of a highly ordered hydrate structure of the polysaccharide, causing a large entropic stress despite the stability brought by stabilization of the complexes. Heating thus leads to dehydration and release of βGP ligands, accompanied by a transfer of protons from CS to βGP due to the simultaneous decrease in the ionization constant of CS as already suggested by [START_REF] Lavertu | Heat-induced transfer of protons from chitosan to glycerol phosphate produces chitosan precipitation and gelation[END_REF]. By losing their charges and hydrated structures, CS chains eventually tend to self-associate and form a gel. These CS-βGP mixtures appeared clearly as a promising thermo-gelling system for cell injection, remaining liquid at physiological pH and turning rapidly into a gel at body temperature.

However, increasing salt concentrations to reach appropriate values of pH also increases the osmolarity of the mixture up to cytotoxic values [START_REF] Ahmadi | Biocompatibility and gelation of chitosan-glycerol phosphate hydrogels[END_REF]. Hence, strategies were implemented to improve the standard system, mostly by replacing partially or totally βGP by other salts. Using ammonium hydrogen phosphate (AHP, Figure 1), a lower salt concentration was required to reach neutral pHs (7.0 -7.2) resulting in physiological osmolarity (300 ± 30 mOsm.L -1 ) [START_REF] Nair | Development of injectable thermogelling chitosaninorganic phosphate solutions for biomedical applications[END_REF]. Other inorganic phosphate salts without polyol moieties, like dipotassium hydrogen phosphate K2HPO4 [START_REF] Ta | Chitosan-dibasic orthophosphate hydrogel: A potential drug delivery system[END_REF], basic sodium phosphates as NaH2PO4, Na2HPO4, or even Na3PO4 (X. Y. [START_REF] Li | Gel-sol-gel thermo-gelation behavior study of chitosan-inorganic phosphate solutions[END_REF][START_REF] Casettari | Characterization of the interaction between chitosan and inorganic sodium phosphates by means of rheological and optical microscopy studies[END_REF], have also been studied to produce thermo-sensitive CS formulations. Sodium hydrogen carbonate (SHC) has been studied as a gelling agent, alone [START_REF] Huang | Modification of an injectable chitosan scaffold by blending with NaHCO 3 to improve cytocompatibility[END_REF][START_REF] Liu | Smart gelation of chitosan solution in the presence of NaHCO3for injectable drug delivery system[END_REF] or in combination with βGP [START_REF] Assaad | Injectable thermosensitive chitosan hydrogels with controlled gelation kinetics and enhanced mechanical resistance[END_REF][START_REF] Ceccaldi | Optimization of Injectable Thermosensitive Scaffolds with Enhanced Mechanical Properties for Cell Therapy[END_REF][START_REF] Deng | Enhanced gelation of chitosan/β-sodium glycerophosphate thermosensitive hydrogel with sodium bicarbonate and biocompatibility evaluated[END_REF]Alinejad, Adoungotchodo, Hui, Zehtabi, & Lerouge, 2018).

Mixtures of CS/AHP [START_REF] Nair | Development of injectable thermogelling chitosaninorganic phosphate solutions for biomedical applications[END_REF] or CS/βGP/SHC [START_REF] Ceccaldi | Optimization of Injectable Thermosensitive Scaffolds with Enhanced Mechanical Properties for Cell Therapy[END_REF]Alinejad, Adoungotchodo, Hui, Zehtabi, & Lerouge, 2018) have been evaluated in vitro as 3D scaffolds for cells. Encapsulation of mouse preosteoblastic cells and fibroblasts as well as human mesenchymal stem cells confirmed the cytocompatibility for days for some formulations.

The feasibility of using these chitosan-based hydrogels as a cell carrier vehicle has therefore been demonstrated. Nevertheless, the numerous variables, pH, osmolarity, salt and chitosan concentrations leading to injectable cytocompatible solutions are all interdependent. A general platform is still lacking, enabling to control independently pH and osmolarity, and therefore gelation kinetics and cytocompatibility, whatever the chitosan concentration. Yet, controlling chitosan concentration is a key to adjust the viscoelastic and microstructural properties of the injected hydrogels, which are of the utmost importance to favor cell anchorage and growth.

In this study we present a rational formulation design to make thermosensitive CS solutions, by partially replacing βGP by AHP in the mixture. With this third component, AHP, we can decouple the control on pH and osmolarity, whatever the CS concentration. We show that these two parameters can be adjusted independently by the simple ratios [salt]/[NH2], where [NH2] denotes the concentration in D-glucosamine monomer in the solution, i.e. the concentration in deacetylated chitin monomers bearing NH2 functions. This approach gives rise to a robust platform, to easily formulate physiological solutions able to form a gel within a couple of minutes at body temperature. We investigated in depth the gelation kinetics for physiological formulations and determined a critical temperature of gelation, under which the system remains liquid overtime, for different CS concentrations. We checked the cytocompatibility of our controlled CS/βGP/AHP formulations and explored the effect of osmolarity on cell activity by encapsulating preosteoblastic mouse cells (MC3T3-E1) and human gingival fibroblasts (HGF).

By means of confocal microscopy and rheology, we showed that hydrogel microstructure and mechanical properties can be varied by changing CS concentration while maintaining fixed pH and osmolarity.

EXPERIMENTAL METHODS

Materials & Cells

Three grades of chitosan (CS) having high purity (impurity level < 1.0%) were purchased from Glentham Life Sciences Limited. All the study was performed with chitosan low molar mass (viscosity 30-100 cps) having the following characteristics given by the supplier: deacetylation degree DD = 91%, as measured by potentiometric titration, average molar mass Mw = 250 kg.mol -1 , as measured by viscosimetry. For pH titration experiments, two other grades with DD of 95 % were also used: chitosan very low molar mass (10 cps) with Mw = 30 kg.mol -1 , and chitosan (100-300 cps) with Mw = 890 kg.mol -1 .

β-glycerophosphate disodium salt hydrate BioUltra (βGP) and ammonium hydrogen phosphate BioUltra (AHP) were purchased from Sigma-Aldrich. Glacial acetic acid (AcOH, 99%) was purchased from Alfa Aesar. Alexa Fluor 488 was purchased from Abcam. The Live/Dead TM kit and AlamarBlue TM reagent were purchased from Invitrogen. MC3T3-E1 pre-osteoblastic cells were procured from ATCC ® , and human gingival fibroblasts were isolated from healthy donors and kindly donated by Laboratory BioS, Reims University. α-ΜΕΜ, DMEM, phosphate buffer saline without calcium and magnesium (PBS) were purchased from Gibco. Deionized distilled water from a Milli-Q system was used to prepare all aqueous solutions.

Methods

2.1.Preparation and sterilization of reactants

Aqueous suspensions of CS were sterilized in autoclave at 121°C for 10 min using standard liquid cycle according to San [START_REF] San Juan | Degradation of chitosan-based materials after different sterilization treatments[END_REF]. After autoclaving, acetic acid (AcOH) was added to the suspension with molar ratio [AcOH]/[NH2] = 1, which was stirred at room temperature (21±1°C) for 24 h. The concentration of the CS stock solution was 1.33 wt%. CS stock solutions were stored at 5°C until further use. βGP, AHP and AcOH solutions at 1 mol.L -1 were prepared in deionized water and sterilized by filtration using syringe filters having pore size of 0.2 µm.

pH titration of CS/βGP solution with AHP solution

Stock solutions of βGP and AHP as well as CS solution at 0.8 wt% were chilled in an ice bath prior to mixture. All pH titrations were performed at 4 ± 1°C. The βGP solution was added dropwise to the CS solution under vigorous stirring. Solutions with molar ratio [βGP]/[NH2] ranging from 0 to 4 were prepared. For each [βGP]/[NH2] ratio, titration was carried out with AHP solution until the pH reached 7.4.

Osmolarity measurement

Osmolality of different solutions was measured with a cryoscopic osmometer (Osmomat 30, Gonotec). The apparatus measures the osmolality (mOsm.kg -1 ), which is proportional through the density to the osmolarity (mOsm.L -1 ). In this study, since all solutions are diluted, the difference between osmolarity and osmolality is insignificant.

Gelation time measurement

The gelation times were determined by the test tube inversion method. The sol phase was defined as the flowing phase and the gel phase as the non-flowing phase when the test tube was inverted.

1 mL of CS/βGP solutions was introduced into 1.5 mL tubes that were incubated in a thermostat at fixed temperature (15-37°C). Time between inversions was as follows: 30 s from 0 to 10 min, 1 min from 10 to 30 min, 5 min from 30 to 60 min, 10 min from 60 to 120 min, 15 min from 2 h to 5 h, 30 min from 5 h to 10 h and 3 h after 10 h. The time at which the gel did not flow anymore was recorded as the gelation time and the time interval between two inversions gave the experimental uncertainty. Each measure was duplicated. The shortest measurable gelation time is given by the duration of thermal equilibration, which in these tubes was measured to be 1 min for a solution at 4°C immersed in a thermostat at 37°C.

Microstructure study of CS/βGP/AHP hydrogels

CS was labelled with Alexa Fluor 488 (Alexa Fluor monomer units:NH2 = 1:150) following a protocol adapted from [START_REF] Zhang | Arginine conjugation affects the endocytic pathways of chitosan/DNA nanoparticles[END_REF]. The fluorescently labelled CS was then added to regular CS suspensions before the autoclave process. Labelled CS chains represented only 0.5 wt% of the total amount of CS in each formulation. The amount of grafted fluorophore was three orders of magnitude smaller than in Zhang et al. (AlexaFluor/NH2 units 1:30000). Hence, the perturbation caused by the labelling was considered to be negligible. Accordingly, there was no difference in gelation time at 37°C between regular CS and fluorescently labelled CS.

A sample holder was filled with a few hundreds of microliters of fluorescently labelled CS/βGP/AHP solution, sealed to prevent evaporation and incubated in a regular oven at 37°C for 30 min. Hydrogels at CS concentrations between 0.5 and 0.8 wt% were observed at room temperature using Laser Confocal Scanning Microscopy (LCSM, Zeiss Axiovert 200M). For each sample, stacks of 10 images of (225 µm x 225 µm) with an interval of 1 µm were recorded in eight random positions at 100 µm depth.

Deconvolution of the image stacks was performed using the Iterative Deconvolve 3D plugin with a theoretical point spread function PSF calculated by the Diffraction PSF 3D plugin (ImageJ).

Deconvoluted images were binarized using Li's method (C. H. [START_REF] Li | An iterative algorithm for minimum cross entropy thresholding[END_REF]. The polymer-rich zones are presented in white and the polymer-poor zones (pores) are presented in black. The porosity was determined as the percentage of black pixels over the total number of pixels. The pore size distribution was characterized quantitatively using a customized algorithm (Python) based on the morphological sieve technique [START_REF] Wu | Pore size distribution in tablets measured with a morphological sieve[END_REF]. This technique calculates on each image the fraction of porosity that is able to fit a structuring element having a given shape and size, here a disk having a given diameter. For each formulation, the average porosity and average pore size distribution were calculated over 8 stacks of 10 images.

Rheological characterization

Rheological properties of physiological formulations prepared with 0.5 wt% or 0.8 wt% of CS were studied using rheometer MCR502 (Anton Paar) with a cone-and-plate geometry (CP50-2TG). The temperature was managed by a Peltier controller and the geometry was previously stabilized at 18°C. Formulations were loaded immediately after preparation and silicone oil was added to the interface between the sample and the air to prevent evaporation.

For isothermal measurement, the temperature was rapidly increased 37°C at a rate of 20°C.min -1 .

The values of storage and loss moduli (G', G") were measured during 6 h, at a constant oscillatory strain of 0.1% and an angular frequency of 1 rad.s -1 . For temperature sweep measurement, the temperature was increased from 4°C to 60°C then decreased from 60°C to 4°C, at a rate of 5°C/min. The values of storage and loss moduli (G', G") were measured during one heating/cooling cycle, at a constant oscillatory strain of 0.1% and an angular frequency of 1 rad.s - 1 .

Injectability evaluation

Injectability test was performed using the protocol reported by [START_REF] Kondziolka | Injection parameters affect cell viability and implant volumes in automated cell delivery for the brain[END_REF]. For visualization, a small amount of rose Bengal dye was added to CS/βGP/AHP mixture. Immediately following the preparation, the mixture was loaded into a 5 mL syringe equipped with a 21G needle and placed on a pump syringe. A volume of 1.5 mL was then injected at a rate of 250 µL.min -1 into a PBS solution at pH 7.4 and 37°C.

Cell encapsulation in thermosensitive chitosan hydrogels

We compared two formulations of hydrogels for cell encapsulation: CS/βGP (600 mOsm.L 1 ) and CS/βGP/AHP (300 mOsm.L -1 ) at fixed CS concentration (0.8 wt%) and at pH 7.4 on two cell models: primary human gingival fibroblast (HGF) and pre-osteoblast immortalized murine cell line (MC3T3-E1). HGF were cultured in DMEM and MC3T3-E1 cells were cultured in α-ΜΕΜ, both cell culture medium supplemented with 10%v/v fetal calf serum and 1%v/v penicillin-streptomycin. For both cell types, the encapsulation protocol was as follows: (1) cells were trypsinized and mixed gently with CS/βGP or CS/βGP/AHP solutions at a final concentration of 10 6 cells.mL -1 ; (2) 300 µL of the mixture containing cells were dispensed into a 24-well culture plate (Nunc™ F96 MicroWell™) and incubated for 1 h at 37°C and 5% CO2 to induce gelation; (3) 1 mL of culture medium supplemented with 10%v/v fetal calf serum and 1%v/v penicillin-streptomycin was added into each well. Hydrogels without cells were also prepared by replacing the cell suspension by culture medium.

Live/Dead assay

After 1 h or 24 h, the encapsulated cells were stained using calcein AM/ethidium homodimer (CalAM/EthD-1) mixture at optimized concentrations ([CalAM]/[EthD-1]: 2 µM/2 µM and 1 µM/0.5 µM for MC3T3-E1 cells and HGF, respectively). For positive control, cells were seeded directly into the well at a density of 10 4 cells.cm -2 . For negative control, cells were killed using methanol solution at 70 wt%. 1 mL of the staining mixture was added into each well and the plates were incubated at 37°C and 5% CO2 for 45 min then the mixture was removed and the hydrogels or cells were rinsed with PBS. 500 µL of PBS solution were added into each well to prevent drying. The distribution and the state of the cells were visualized using a laser scanning confocal microscope (LSM 710, Zeiss) using lasers at 488 nm and 561 nm. Each hydrogel sample was scanned from bottom to 300 µm depth at 8 random positions. Stacks were taken at 20x magnification and a 512x512 pixel resolution with the pinhole at 1 Airy unit and interval of 4 µm. Recorded images were then analyzed with ImageJ using 3D Object counter plugin [START_REF] Bolte | A guided tour into subcellular colocalization analysis in light microscopy[END_REF]. The cell viability (%) was determined as the number of living cells over the total number of cells. Experiments were carried out with n = 2 for MC3T3 and n = 1 for HGF and all samples were prepared in triplicate.

AlamarBlue assay

The cell metabolism was evaluated using the AlamarBlue assay at 1 h and 24 h following encapsulation. After incubation, a solution of 10 %v/v AlamarBlue in culture medium (1 mL/well) was added directly on each hydrogel and incubated for 3 h at 37°C and 5 % of CO2.

Then, the supernatant was homogenized by pipetting and 150 µL were transferred into a 96-well plate (Nunc™ F96 MicroWell™). The intensity of fluorescence was determined by using Xenius XM spectrofluorometer (SAFAS Monaco) at 560 nm excitation and 590 nm emission. A value of the metabolic activity per cell was estimated by normalizing the measured intensity by the number of living cells obtained from the Live/Dead assay.

RESULTS & DISCUSSION

In all the following, [NH2] denotes the concentration in D-glucosamine monomer in the solution.

Controlling the pH using CS/βGP/AHP formulations

Control of the pH of thermogelling CS/βGP solutions is determinant for the viability and the injectability of encapsulated cells. Cell viability requires that the pH be maintained close to 7.4 during cell suspension and storage until injection and immersion in other buffering media. The gelation time of CS/βGP solutions was also shown to depend greatly on pH [START_REF] Chenite | Rheological characterisation of thermogelling chitosan/glycerol-phosphate solutions[END_REF].

Accordingly, we measured that the gelation time of CS/βGP solutions increases from 1 min to 10 h when pH decreases from 7.5 to 7.1 (Figure S1). In the present study, CS was dissolved in acidic aqueous solutions with a fixed molar ratio [AcOH]/[NH2] = 1. This ratio was chosen as the minimum amount of acid required to ensure the solubilization of chitosan chains [START_REF] Rinaudo | Chitin and chitosan: Properties and applications[END_REF]. In CS/βGP/AHP formulations, βGP and AHP are both weak bases that raise the pH of the solution. As demonstrated thereafter, the pH can be finely controlled by the molar ratios between Those iso-pH curves are linear in the studied range of pH from 7.0 to 7.4. Formulations using the ratios of the iso-pH curve at 7.4 but different initial CS concentrations (0.5 and 0.6 wt%) exhibit a pH 7.4 as well, which confirms that the pH is fully determined by the ratios [βGP]/[NH2] and

[AHP]/[NH2], irrespectively of the CS concentrations.

The pH of CS and CS/βGP solutions is known to decrease substantially with increasing temperature [START_REF] Lavertu | Heat-induced transfer of protons from chitosan to glycerol phosphate produces chitosan precipitation and gelation[END_REF][START_REF] Supper | Rheological study of chitosan/polyol-phosphate systems: Influence of the polyol part on the thermo-induced gelation mechanism[END_REF]. A similar effect was observed in the CS/βGP/AHP formulations. For instance, the pH of a solution at CS 0.8 wt% decreases from 7.40 ± 0.05 at 4°C to 6.95 ± 0.05 at 37°C (Figure S2). Therefore, it seems more appropriate to target formulations having a pH 7.4 at 4°C to ensure that the pH of cell-laden solutions remains within a physiological range (typically 7.0-7.4) during cell dispersion and encapsulation.

NH ⇌ NH H , KCS (1)
AcOH ⇌ AcO H , KAc (2) AHP ⇌ AHP H ⇌ AHP 2H , K , & K , ( 3
) βGP ⇌ βGP H ⇌ βGP 2H , K , & K , ( 4
) H O ⇌ H OH (5)
The values of the dissociation constants of AHP and βGP were measured by pH titration at 4°C and 25°C (Figure S3). Only the second equilibrium is observed for pH ranging from 2 to 8: K , = 10 -6.7 and K , = 10 -6.2 . For CS and AcOH, values are chosen from the literature: KCS = 10 -7.1 (4°C) (Filion et al., 2013) and KAc = 10 -4.8 [START_REF] Kortum | Dissociation Constants of Organic Acids in Aqueous Solution[END_REF].

To a first approximation, we assume here that CS behaves like a small molecule and we do not take into account its polyelectrolyte nature. Writing electroneutrality, conservation of matter and acid-base equilibria, we find that for the pH range of interest (7.0-7.4), the ratios The model predicts the experimental values of the slopes with a good precision including the slight increase when the pH decreases from 7.4 to 7.0, as shown in Figure 2b and Table S1.

However, it fails to predict the intercept. We attribute this discrepancy to the polyelectrolyte nature of CS. In particular, the protonation of NH2 along the CS chains is very sensitive to entropic effects and pKCS increases with increasing salt concentration (Filion et al., 2007).

Accordingly, predictions of pH for CS/βGP accounting for these polyelectrolyte effects using a serve as an efficient platform to prepare CS/βGP/AHP formulations having a controlled pH.

CS/βGP/AHP formulations with adjustable pH and osmolarity

The contribution of βGP, AHP and AcOH to the osmolarity of the formulations was determined by measuring independently their osmolarity at different concentrations. We verified that osmolarity was proportional to the molar concentration for all three components (Figure S5). As expected from the dissociation constants, the osmolarity induced by βGP is higher than the one induced by AHP. The osmolarity of a CS solution at 0.8 wt% is negligible. The total osmolarity (expressed in Osm.L -1 ) was calculated as the linear combination of the three contributions:

Osmolarity = 1.03⋅[AcOH] + 2.35⋅[βGP] + 1.82⋅[AHP] (7)
Normalizing Equation ( 7) by [NH2] gives a relationship between [AHP]/[NH2] and [βGP]/[NH2]:

20 [ ] [$ % ] = -A . [ ] [$ % ] -A . [ 01 ] [$ % ]
A .

1:;<=>?@AB [$ % ]

(8)

where A1 = 1.29, A2 = 0.57 and A3 = 0.55 are dimensionless constants.

For a fixed ratio [AcOH]/[NH2], we obtain linear iso-osmolarity lines between [AHP]/[NH2] and

[βGP]/[NH2]. Unlike the iso-pH curves above, the intercept of these curves depends on the CS concentration. The curves obtained for a physiological osmolarity of 300 mOsm.L -1 are presented in Figure 3a for CS concentrations ranging from 0.4 to 1.2 wt%. For a given CS concentration, the intersection between the iso-pH (pH = 7.4) and iso-osmolarity lines give the unique formulation satisfying both conditions on targeted pH and osmolarity.

The representation in Figure 3a provides a way to determine the formulations of interest for thermosensitive cell-laden hydrogels. In particular, it shows that the double constraint on pH and osmolarity restricts the range of pertinent CS concentrations. As shown in Figure 3a, CS concentration of 0.4 wt% is the lowest concentration allowing to achieve physiological pH (7.4) and osmolarity (300 mOsm.L -1 ). CS concentrations lower than 0.4 wt% would require either a hypotonic mixture to satisfy pH = 7.4 or a basic pH (>7.4) to satisfy an osmolarity of 300 mOsm.L -1 . We verified that the physiological formulation at 0.4 wt% corresponding to

[βGP]/[NH2] = 5.5 and [ΑΗP]/[NH2] = 0 forms a macroscopically homogeneous gel after 2 min at 37°C (Figure 3b). In Figure 3a, the highest CS concentration allowing to achieve physiological pH and osmolarity is 1.2 wt%. However, this concentration (corresponding to

[βGP]/[NH2] = 0) is not pertinent for practical applications, because the gelation time at room temperature (i.e. before injection) decreases when decreasing [βGP]/[NH2], as will be presented in section 3. We found that a minimum ratio [βGP]/[NH2] = 1.0, corresponding to a CS concentration of 0.9 wt%, is required to maintain the thermosensitive solution in a liquid state for out cell trypsination, dispersion in CS/salt solution and injection. Accordingly, solutions prepared at a higher concentration than 0.9 wt% form a gel in less than 1 h at room temperature, as shown These results and methodology show how the pH and osmolarity of CS/βGP/AHP mixtures can be independently adjusted for different CS concentrations. It is worth to note that the concentration range for cytocompatible values of pH and osmolarity has been measured here for one type of CS having a given degree of deacetylation (DD = 91 %) and a given molar mass distribution (viscosity 30-100 cps, Mw = 250 kg.mol -1 ). It has been reported that the gelation rate of CS/βGP solutions at given pH increases with CS of higher DD (Ruel-Gariépy, [START_REF] Ruel-Gariépy | Characterization of thermosensitive chitosan gels for the sustained delivery of drugs[END_REF]. Therefore, the upper bound of the accessible concentration range is expected to decrease with increasing DD. As regards the molar mass, we verified that the pH can be adjusted by the salts/NH2 ratios whatever the chitosan chain length is as long as the deacetylation degree is known, as shown for CS with DD = 95 % for two different molar masses (Figure 4). Hence, this rationalized formulation is applicable to other types of CS.

Gelation kinetics of CS/βGP/AHP formulations

It is important to notice that the gelation time of CS/βGP systems depends greatly on the pH [START_REF] Chenite | Rheological characterisation of thermogelling chitosan/glycerol-phosphate solutions[END_REF]. As a result, we focus here on the gelation kinetics at a fixed pH = 7.4.

Figure 5 shows the evolution of storage, G', and loss modulus, G", during a heating-cooling cycle between 4 and 60°C for CS/βGP and CS/βGP/AHP, both at pH 7.4. While heating, both moduli exhibit a sharp increase and the gel point (G' = G") is reached near 45°C for CS/βGP and 40°C

for CS/βGP/AHP. Upon cooling, a large hysteresis is observed and both moduli decrease slowly below 30°C. This weak thermoreversibility is very similar to what has been reported earlier by [START_REF] Chenite | Rheological characterisation of thermogelling chitosan/glycerol-phosphate solutions[END_REF]. A more systematic characterization of the thermosensitivity of CS/βGP/AHP mixtures was performed using the test tube inversion method (Figure S6a). We verified that the time measured by this method is consistent with the time to reach the gel point as determined by rheology (Figure S6b-c). In a first series of experiments, we studied the effect of the relative salt ratio

[AHP]/[βGP] on the gelation time. For that, a fixed CS concentration of 0.8 wt% was chosen.

Four ratios [βGP]/[NH2] were studied (0, 1.5, 3 and 5.5) and the corresponding [AHP]/[NH2] ratios were adjusted to maintain pH = 7.4 (1.9, 1.4, 0.8 and 0, respectively). In those solutions, the osmolarity varied from 300 to 600 mOsm.L -1 . Figure 6a shows that the gelation time decreases by at least four orders of magnitude when temperature increases from 15°C to 37°C for all four formulations. At a given temperature, the gelation time decreases with decreasing Similar experiments were performed in a second series with formulations defined by the iso-pH (7.4) and iso-osmolarity (300 mOsm.L -1 ) curves for different CS concentrations ranging from 0.5 to 0.8 wt%. The gelation time decayed strongly with increasing temperature in a similar way to the first series, as shown in to maintain physiological conditions). Again, a critical power law behavior was observed as shown in Figure 6e (and Figure S7c-d). The critical temperature, and therefore the thermal stability of the solution, decreased with increasing CS concentration.

From a practical point of view, these results show that all the formulations having physiological pH and osmolarity present an adequate gelation kinetics for use as cell delivery vehicles. At 21°C, they remain liquid for more than 2 h, which is appropriate for cell dispersion and manipulation. At body temperature, they form a gel within less than 1 min, therefore preventing cells from being dispersed away. The injectability of the formulation having a CS concentration of 0.8 wt% was verified using standard injection parameters into a PBS solution at 37°C. A continuous injection could be performed with the quasi-instantaneous formation of a macroscopic hydrogel (Video S1).

Cell encapsulation in CS/βGP/AHP hydrogels

The viability and the metabolism of encapsulated cells was studied using Live/Dead and AlamarBlue assays with formulations having a CS concentration of 0.8 wt% and pH 7. ) and primary human gingival fibroblasts (HGF). Cells were dispersed in the liquid formulations at 20°C at a concentration of 10 6 cells. mL -1 . Encapsulation was achieved by gelation at 37°C for 1 h in an incubator prior to the addition of culture medium. Cells were stained and observed 1 h and 24 h after encapsulation. The staining efficiency was verified on positive and negative controls (Figure S8). The homogeneity of the cell dispersion was verified by cell counting in multiple areas of the gels (Figure S9) and a local cell concentration of 1.2 x 10 6 ± 0.2 cells per cm 3 was found, which is in agreement with the density of seeded cells.

Typical images obtained with Live/Dead staining are shown in Figure 7a-b for MC3T3-E1 and HGF, respectively. At 1 h, the viability for both cell types is high (90%) for the two studied formulations (Figure 7c). At 24 h, differences are observed: for MC3T3-E1 cells, a significantly better cell viability is observed in CS/βGP/AHP hydrogel (78 ± 6 %) than in CS/βGP hydrogel (52 ± 5 %); for HGF, cell viability remains high in both hydrogels. The metabolic activity of the encapsulated cells was assessed with an AlamarBlue assay (Figure 7d). A value of the metabolic activity per cell was estimated by normalizing the measured intensity by the number of living cells. At 1 h, a clear difference can be seen between the two cell types: MC3T3 have a much higher metabolic activity in the isotonic formulation CS/βGP/AHP than in CS/βGP while HGF have a similar metabolic activity in both formulations. At 24 h, the same trend is observed but the metabolic activity drops significantly for both cell types. These differences in viability and activity may be explained by the different ability of the cell types to adapt to hypertonic conditions [START_REF] Dmitrieva | Hypertonic stress response[END_REF]. To our knowledge, no study has been reported on MC3T3-E1 whereas it has been shown that HGF are resistant to osmotic shocks [START_REF] Huynh | Rinsing with saline promotes human gingival fibroblast wound healing in vitro[END_REF]. In our study, the encapsulated cells exhibit rounded phenotypes after 24 h, which reflects a poor adhesion to the chitosan matrix. This observation is in agreement with a previous study on encapsulated MC3T3-E1 in CS/AHP hydrogels [START_REF] Nair | Development of injectable thermogelling chitosaninorganic phosphate solutions for biomedical applications[END_REF]. Therefore, the proposed approach could be advantageously combined with bio-functionalization strategies providing adhesion sites to the encapsulated cells. For example, covalent grafting of peptides or proteins responsible for cell adhesion on chitosan chains have shown promising results [START_REF] Custodio | Immobilization of fibronectin in chitosan substrates improves cell adhesion and proliferation[END_REF]Tsai, Chen, & Liu, 2013).

Microstructure and mechanics of CS/βGP/AHP hydrogels

The possibility to vary the CS concentration while maintaining physiological conditions, provides a way to adjust the microstructure and the mechanical properties of cytocompatible formulations.

The microstructure of hydrogels (pH = 7.4 and 300 mOsm.L -1 ) was characterized by LSCM using fluorescently labelled formulations. After 30 min at 37 °C, heterogeneous structures were observed composed of interconnected polymer-rich and polymer-poor phases, the latter being qualified as pores (a). We measured that the porosity (%) decreases linearly with increasing CS concentration as shown in Figure 8b, from 75 ± 5 % for CS 0.5 wt% to 63 ± 7 % for CS 0.8 wt%.

Such a linear relationship was reported earlier for CS/βGP [START_REF] Crompton | Morphology and gelation of thermosensitive chitosan hydrogels[END_REF]. Besides porosity, the pore size is a crucial parameter governing cell-cell interactions and vascularization [START_REF] Kaivosoja | Chemical and physical properties of regenerative medicine materials controlling stem cell fate[END_REF]. Targeted pore sizes vary from tens to hundreds of micrometers depending on the cell type and on the application [START_REF] Kaivosoja | Chemical and physical properties of regenerative medicine materials controlling stem cell fate[END_REF][START_REF] O'brien | The effect of pore size on cell adhesion in collagen-GAG scaffolds[END_REF]. With the present systems, we observed large interconnected elongated macropores. The longest dimension of the largest pores was found to reach 15 µm for CS 0.8 wt% and up to 30 µm for CS 0.5 wt%. Pore size distribution was analyzed quantitatively by image analysis using a morphological sieve algorithm. The distributions in Figure 8c show the fraction of porosity able to fit a circular structuring element of a given diameter. This analysis the ones of most biological tissues. For this reason, these hydrogels might be more particularly appropriate for very soft tissues, such as brain tissue which G' is in the range 100-300 Pa [START_REF] Tabet | Mechanical Characterization of Human Brain Tissue and Soft Dynamic Gels Exhibiting Electromechanical Neuro-Mimicry[END_REF]. Accordingly, CS gels were shown to be suitable matrices for cells of the nervous system [START_REF] Kornev | Hydrogel-assisted neuroregeneration approaches towards brain injury therapy: A state-of-the-art review[END_REF].

CONCLUSION

In this study, we report a rational formulation design to prepare thermo-gelling systems composed of chitosan CS and two phosphate salts, βGP and AHP. We show how the pH and the osmolarity can be adjusted independently by the ratios [salts]/[NH2], arising possibility to obtain physiological formulations having tunable CS concentration. Moreover, a systematic study of the gelation kinetics revealed that the gelation time follows a power law that diverges at a critical temperature. With this approach, injectable solutions having a pH of 7.4 and a physiological osmolarity of 300 mOsm.L -1 can be produced for CS concentrations ranging from 0.4 to 0.9 wt%.

The cytocompability of the subsequent gels was demonstrated by in vitro encapsulation of MC3T3-E1 and HGF cells. Elastic moduli of the gels can reach up to 300 Pa, a value found in brain tissues. The proposed methodology could be applied to solutions of chitosan or functionalized chitosan having various molar mass and deacetylation degrees, therefore providing a general platform to prepare valuable thermosensitive hydrogels for cell encapsulation and delivery.

acknowledged. 547

Figure 1 .

 1 Figure 1. Structures of chitosan CS, β-glycerophosphate βGP and ammonium hydrogen phosphate AHP.

  each phosphate salt and D-glucosamine monomer, [βGP]/[NH2] and [AHP]/[NH2]. We carried out pH titration at 4°C by adding AHP (1 mol.L -1 ) to solutions having an initial CS concentration of 0.8 wt% and different ratios of [βGP]/[NH2]. The evolution of pH as a function of [AHP]/[NH2] is reported in Figure 2a. For a given [βGP]/[NH2], all titration curves show an increase in pH with increasing [AHP]/[NH2]. Similarly, for a given [AHP]/[NH2], pH increases with increasing [βGP]/[NH2]. From those curves, we extracted the pairs of ratios [βGP]/[NH2] and [AHP]/[NH2] producing a given pH and built the iso-pH curves represented in Figure 2b.

Figure 2 .

 2 Figure 2. (a) Evolution of the pH as a function of [AHP]/[NH2] for formulations with initial CS concentration of 0.8 wt% and [βGP]/[NH2] ranging from 0 to 4. (b) Iso-pH curves at 7.0, 7.2 and 7.4 as a function of [AHP]/[NH2] and [βGP]/[NH2]. Experimental data are shown as filled symbols. Linear fitting curves are shown as continuous lines and model predictions as interrupted lines. Predictions accounting for the polyelectrolyte effect based on Filion et al., 2007 are shown as open symbols for CS/βGP mixtures ([AHP]/[NH2] = 0).

  [AHP]/[NH2] and [βGP]/[NH2] satisfy the following relation (see detailed calculation in SI): ], [βGP], [AHP], [NH ] are respectively the initial concentrations in AcOH, βGP, AHP and D-glucosamine monomer, DD is the degree of deacetylation of CS. Factors CAc, CβGP, CAHP and CCS are all functions of the concentration in protons h as follows: ratio [AcOH]/[NH2], here equal to 1, equation (6) captures that for a given pH, the relationship between [AHP]/[NH2] and [βGP]/[NH2] is linear. The negative slope -( )*+ ( ,-+ depends weakly on the pH of the mixture and the dissociation constants of the two salts βGP and AHP.

  model by Filion et al. (Figure S4) are in perfect agreement with our experimental measurements, as indicated by the open symbols in Figure 2b. The simple small molecule model also confirms that the linear relationship between [AHP]/[NH2] and [βGP]/[NH2] is independent on the CS concentration as long as the ratio [AcOH]/[NH2] is fixed. Therefore, such iso-pH curves may

  for example for a 1 wt% solution ([βGP]/[NH2] = 0.7, [ΑΗP]/[NH2] = 1.7) which formed a gel within 30 min (Figure3c).

Figure 3 .

 3 Figure 3. (a) Iso-osmolarity curves at 300 mOsm.L -1 (black) and iso-pH curve at 7.4 (blue) as a function of [AHP]/[NH2] and [βGP]/[NH2] for CS concentration ranging from 0.4 to 1.2 wt%. Shaded area indicates formulations for which gelation time at 21°C is less than 1 h. (b) Hydrogel (300 mOsm.L -1 ; pH = 7.4) obtained with 0.4 wt% of CS after 2 min of incubation at 37°C, (c) Hydrogel (300 mOsm.L -1 ; pH = 7.4) obtained with 1.0 wt% of CS after 30 min at 21±1°C.

Figure 4 .

 4 Figure 4. Evolution of the pH as a function of [AHP]/[NH2] in formulations with initial CS concentration of 0.8 wt% and fixed [βGP]/[NH2] = 1.5, 2.5 and 3.5 for three different chitosan grades (10 cps, 30-100 cps and 100-300 cps) having DD = 95%, 91% and 95%, respectively.

Figure 5 .

 5 Figure 5. Evolution of the storage and loss modulus during a heating/cooling cycle for CS/βGP mixture (a) and CS/βGP/AHP mixture (b), both at fixed pH of 7.4, CS concentration of 0.8 wt%.

[

  Figure5supports the hypothesis that the polyol part of βGP hinders the interactions between CS chains as suggested by[START_REF] Supper | Rheological study of chitosan/polyol-phosphate systems: Influence of the polyol part on the thermo-induced gelation mechanism[END_REF].

Figure 6d .

 6d For a given temperature, gelation occurred faster with increasing CS concentration (together with decreasing [βGP]/[NH2] and increasing [AHP]/[NH2]

Figure 6 .

 6 Figure 6. (a-c) Gelation kinetics of formulations having fixed CS 0.8 wt% and pH 7.4 and various [βGP]/[NH2] and [AHP]/[NH2] ratios. Gelation time as a function of temperature (a) and (T-Tc) (b). Full lines show fits by a power law. (c) Critical temperature as a function of [βGP]/[NH2] and [AHP]/[NH2]. (d-f) Gelation kinetics of physiological formulations (pH = 7.4 and 300 mOsm.L -1 ) for CS concentration ranging from 0.5 to 0.8 wt%. Gelation time as a function of temperature (d) and (T-Tc) (e). Full lines show fits by a power law. (f) Critical temperature as a function of CS concentration.

  4. Two formulations were compared to assess the effect of osmolarity: an isotonic formulation denoted CS/βGP/AHP (300 mOsm.L 1 ) produced with ratios [βGP]/[NH2] = 1.5 and [AHP]/[NH2] = 1.4, and a hypertonic formulation denoted CS/βGP (600 mOsm.L -1 ) obtained with ratios [βGP]/[NH2] = 5.5 and [AHP]/[NH2] = 0. Two cell types were studied: pre-osteoblast immortalized murine cell line (MC3T3-E1

Figure 7 .

 7 Figure 7. Cell encapsulation in CS/βGP/AHP and CS/βGP hydrogels for 1 h and 24 h. (a-b) Confocal images of Live/Dead assays for MC3T3 (a) and HGF (b) cells. Images are obtained from the projection of 76 slices taken at a random position in the hydrogel with a field of 425 x 425 x 300 µm. (c) Cell viability. (d) Normalized metabolic activity. Statistically significant at p < 0.05 (*) or p < 0.0001 (****).

Figure 9 .

 9 Figure 9. Storage modulus (G') and loss modulus (G") as function of time at 37°C for physiological formulations prepared with a CS concentration of 0.5 wt% ([βGP]/[NH2] = 4.0 and [AHP]/[NH2] = 0.5) and 0.8 wt% ([βGP]/[NH2] = 1.5 and [AHP]/[NH2] = 1.4).

2.11. Statistical analysisStatistical analysis was performed using a t-test with normality assumption. Statistical significance was accepted at the level of p < 0.05.

ACKNOWLEDGMENTS

We thank Louison Blivet for technical support on the study of Mw effect, Michel Cloître and Mickael Pomes-Hadda (C3M, ESPCI Paris-PSL) for their help and advice on rheology, Jean Baudry and Florence Condamine (CBI, ESPCI Paris-PSL) for their help on osmometry, Lamia El Guermah and Remy Agniel (ERRMECe, CYU) for technical support on the cell culture and characterization. We thank Etienne Decencière (CMM, MINES ParisTech-PSL) for his help and advice on image analysis. We thank the Laboratory BioS, Reims University for their donation of

clearly reveals that the characteristic size of the macropores increases with decreasing CS concentration. Indeed, 98% of the porosity can contain structuring elements having diameters ranging from 1 to 9 µm for CS 0.8 wt% and from 1 to 14 µm for CS 0.5 wt%. The mechanical properties of hydrogels prepared with 0.5 and 0.8 wt% at 4°C were characterized by oscillatory rheology during an isotherm at 37°C as shown in Figure 9. Gelation occurs within one minute for both formulations with a rapid increase in the storage modulus (G') reaching approximately 30 Pa for 0.5 wt% and 180 Pa for 0.8 wt% at 1 h of incubation. After 6 h, G'


Pascal Ochem
email: ochem@lirmm.fr

A family of formulas with reversal of arbitrarily high avoidability index

We present a family of avoidable formulas with reversal whose avoidability index is unbounded. We also complete the determination of the avoidability index of the formulas with reversal in the 3-avoidance basis.

Introduction

The notion of formula with reversal [START_REF] Currie | A family of formulas with reversal of high avoidability index[END_REF][START_REF] Currie | Avoidance bases for formulas with reversal[END_REF] is an extension of the notion of classical formula such that a variable x can appear both as x and x R with the convention that in an occurrence h of the formula, h(x R ) is the reverse (i.e., mirror image) of h(x). For example, the word w = 20012210101122201 contains the occurrence h : x → 01, y → 221 of the formula xyx • xy R • x R because h(xyx) = 0122101, h(xy R ) = 01122, h(x R ) = 10 are all factors of w. The avoidability index λ(F ) of a formula with reversal F is the minimum number of letters contained in an infinite word avoiding F .

Currie, Mol, and Rampersad [START_REF] Currie | A family of formulas with reversal of high avoidability index[END_REF] have asked if there exist formulas with reversal with arbitrarily large avoidability index. They considered the formula

ψ k = xy 1 y 2 • • • y k x • y R 1 • y R 2 • • • y R k and obtained that λ(ψ 1 ) = 4, λ(ψ 2 ) = λ(ψ 3 ) = λ(ψ 6 ) = 5, 5 λ(ψ 4 ) 6, 5 λ(ψ 5 ) 7, 4 λ(ψ k ) 6 if k 7
and k ≡ 0 (mod 3), and 4 λ(ψ k ) 5 if k 9 and k ≡ 0 (mod 3). They conjecture that λ(ψ k ) = 5 for all k 2. Computational experiments suggest that the upper bound λ(ψ k ) 5 for k 3 is witnessed by the image of every 7 4 + -free ternary word under the following (k + 3)-uniform morphism where k = 3t + i, t 1, and 0 i 2.

0 → (012) t+1-i (0123) i 1 → (013) t+1-i (0134) i 2 → (014) t+1-i (0142) i
We do not try to prove that such words actually avoid ψ k . Instead, we give a positive answer to their original question with Theorem 1 below. Consider the formula

φ k = x 0 x 1 •x 1 x 2 •. . .•x k-1 x 0 •x R 0 •x R 1 •. . .•x R k-1 . That is, φ 1 = x 0 x 0 •x R 0 , φ 2 = x 0 x 1 • x 1 x 0 • x R 0 • x R 1 , φ 3 = x 0 x 1 • x 1 x 2 • x 2 x 0 • x R 0 • x R 1 • x R 2 , . . . Theorem 1. For every fixed b, there exists k such that b < λ(φ k ) k + 1.
This result contrasts with the situation of classical formulas (without reversal). Clark [START_REF] Clark | Avoidable formulas in combinatorics on words[END_REF] has shown in 2001 that there exist formulas with index 5 (such as ab.ba.ac.bc.cda.dcd), but no avoidable classical formula with higher index is currently known.

Currie, Mol, and Rampersad [START_REF] Currie | Avoidance bases for formulas with reversal[END_REF] have also determined the 3-avoidance basis for formulas with reversal, which contains the minimally avoidable formulas with reversal on 3 variables. They obtained several bounds on the avoidability index of the formulas with reversal in the 3-avoidance basis. The next two results finish the determination of the avoidability index of these formulas.

Theorem 2. The following formulas are simultaneously 2-avoidable:

• xyzyx • zyxy R z • xyzyx • zy R xyz • xyzyx • zy R xy R z • xyzy R x • zyxy R z • xyzy R x • zy R xyz Theorem 3. The formulas xyzx • yzxy • z R and xyzx • yz R xy are simulta- neously 3-avoidable.
Theorems 1 to 3 are proved in Sections 2 to 4, respectively. A word w is d-directed if for every factor f of w of length d, the word f R is not a factor of w.

Remark 4. If a d-directed word contains an occurrence

h of x • x R , then |h(x)| d -1.
In order to express the simultaneous avoidance of similar formulas, as in Theorems 2 and 3, we introduce the notation x U to represent equality up to mirror image. That is, if h(x) = w, then h(x R ) = w R and h(x U ) ∈ w, w R . For example, avoiding xyxy and xyx R y simultaneously is equivalent to avoiding xyx U y. Notice that the notion of undirected avoidability recently considered by Currie and Mol [START_REF] Currie | The undirected repetition threshold and undirected pattern avoidance[END_REF] corresponds to the case where every occurrence of every variable of the formula is equipped with -U .

Recall that a word is (β + , n)-free if it contains no repetition with exponent strictly greater than β and period at least n. Also, a word is (β + )-free if it is (β + , 1)-free.

The C code to find and check the morphisms in this paper is available at http://www.lirmm.fr/ ~ochem/morphisms/reversal.htm.

Formulas with unbounded avoidability index

Let us first show that for every k 2, φ k is avoided by the periodic word ( 0 1 • • • k ) ω over (k + 1) letters. This word is 2-directed, so every occurrence h of φ k is such that |h(x i )| = 1 for every 0 i < k by Remark 4. Without loss of generality, h(x 0 ) = 0 . This forces h(x 1 ) = 1 , h(x 2 ) = 2 , and so on until h(x k-1 ) = k-1 and h(x 0 ) = k , which contradicts h(x 0 ) = 0 . Thus λ(φ k ) k + 1.

Let b be an integer and let w be an infinite word on at most b letters. Consider the Rauzy graph R of w such that the vertices of R are the letters of w and for every factor uv of length two in w, we put the arc -→ uv in R. So R is a directed graph, possibly with loops (circuits of length 1) and digons (circuits of length 2). Since w is infinite, every vertex of R has out-degree at least 1. So R contains a circuit C i of length i with 1 i b. Let c 0 , c 1 , . . . , c i-1 be the vertices of C i in cyclic order. Let k be the least common multiple of 1, 2, . . . , b. Since i divides k, w contains the occurrence h of φ k such that h(x j ) = c j mod i for every 0 j < k. Thus λ(φ k ) > b.

Formulas that flatten to xyzyx • zyxyz

Notice that avoiding simultaneously the formulas in Theorem 2 is equivalent to avoiding F = xyzy U x • zy U xy U z • y R . The fragment y R is here to exclude the classical formula xyzyx • zyxyz. Indeed, even though Gamard et al. [START_REF] Gamard | Avoidability of circular formulas[END_REF] obtained that λ(xyzyx • zyxyz) = 2, a computer check shows that xyzyx • zyxyz and F cannot be avoided simultaneously over two letters, that is, xyzy U x • zy U xy U z is not 2-avoidable.

We use the method in [START_REF] Ochem | A generator of morphisms for infinite words[END_REF] to show that the image of every

4 . 5 +

 45 The formulas xyzx • yzxy• z R and xyzx • yz R xy Notice that avoiding xyzx • yzxy • z R and xyzx • yz R xy simultaneously is equivalent to avoiding F = xyzx • yz U xy • z R .We use the method in[START_REF] Ochem | A generator of morphisms for infinite words[END_REF] to show that the image of every 7 -free word over Σ 4 under the following 9-uniform morphism is 131 90 + , 28 -free. We also check that such a ternary word is 4-directed.


Géraldine Rios

Caroline Lacoux

Vianney Leclercq

Anna Diamant

Kevin Lebrigand

Adele Lazuka

Emmanuel Soyeux

Sebastien Lacroix

Julien Fassy

Aurélie Couesnon

Kévin Lebrigand

Adèle Lazuka

Sébastien Lacroix

Richard Thiery

Bernard Mari

Christian Pradier

Rainer Waldmann

Pascal Barbry
email: barbry@ipmc.cnrs.fr

Monitoring SARS-CoV-2 variants alterations in Nice neighborhoods by wastewater nanopore sequencing

de documents scientifiques 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 (402 words)

New strategies are needed to track the emergence and spread of SARS-CoV-2 variants.

Epidemiological surveillance has largely relied on individual testing by quantitative real-time polymerase chain reaction (RT-qPCR) or digital droplet polymerase chain reaction (ddPCR).

However, a patient-oriented approach is hardly sustainable for a long period in a large population. Individual testing has also some intrinsic limitations. Information is mostly derived from symptomatic cases and does not necessarily reflect the full dynamics of COVID-19 diffusion. The incomplete identification of asymptomatic, pre-symptomatic or postsymptomatic cases by individual testing can be overcome by analyses of pooled samples that reflect the entire population [START_REF] Pettit | All In': a pragmatic framework for COVID-19 testing and action on a global scale[END_REF]. The principal sites of SARS-CoV-2 replication are the upper and lower respiratory tracts, but the virus also replicates in the intestine, leading to high concentrations of SARS-CoV-2 in human stool [START_REF] Zheng | Viral load dynamics and disease severity in patients infected with SARS-CoV-2 in Zhejiang province, China[END_REF]. In consequence, analysis of SARS-CoV-2 in wastewater appears a promising approach for cost-effective population-scale SARS-CoV-2 epidemiology. Several recent studies that analyzed SARS-CoV-2 titers in Wastewater Treatment Plants (WWTPs) by qPCR revealed a good correlation between SARS-CoV-2 incidence rates and the virus titers in wastewater [START_REF] Medema | Implementation of environmental surveillance for SARS-CoV-2 virus to support public health decisions: Opportunities and challenges[END_REF][START_REF] Peccia | Measurement of SARS-CoV-2 RNA in wastewater tracks community infection dynamics[END_REF][START_REF] Rosa | First detection of SARS-CoV-2 in untreated wastewaters in Italy[END_REF][START_REF] Ahmed | First confirmed detection of SARS-CoV-2 in untreated wastewater in Australia: A proof of concept for the wastewater surveillance of COVID-19 in the community[END_REF][START_REF] Wu | SARS-CoV-2 Titers in Wastewater Are Higher than Expected from Clinically Confirmed Cases[END_REF][START_REF] Wurtzer | Evaluation of lockdown effect on SARS-CoV-2 dynamics through viral genome quantification in waste water[END_REF].

The emergence of variants of concern (VOCs) associated with higher transmissibility and/or immune evasion in England (B.1.1.7), South Africa (B.1.351), South America (P.1), the United States of America (B.1.525) and India (B.1.617.2) recently added a further layer of complexity. Measurement of virus concentrations needs now to be complemented by a direct quantification of the individual VOCs. This can be done with variant-specific RT-qPCR/ddPCR assays, but the continuing emergence of new variants calls for more generic approaches.

High throughput sequencing (HTS) of SARS-CoV-2 overcomes the limitations of qPCR since it allows the identification of both known and novel mutations in the virus sequence. HTS also allows the definition of the relative abundance of individual virus variants in a sample and thus is well suited for the analysis of complex virus populations in samples such as wastewater.

Our aim was to evaluate a wastewater SARS-CoV-2 sequencing approach in a real-life configuration. We performed this study in Nice (France). Samples covering 416,755 inhabitants were collected from the inlet of the wastewater treatment plant and upstream in sewers draining different neighborhoods once a month between October 2020 and March 2021. We show here that nanopore sequencing of wastewater not only revealed the dynamics of known SARS-CoV-2 lineages within the city of Nice but also allowed the identification of novel variants that would have escaped identification with variant-specific qPCR assays.

Results and Discussion (1931 words)

Our goal was to analyze lineage composition in multiple neighborhoods, including specific areas that are more susceptible to outbreaks (e.g. zones with a lower standard of living, higher population density). We collected 113 distinct wastewater samples from 21 different locations spread all over the territory of Nice (Figure 1A-B, Supplementary Table 1) once a month between October 2020 and March 2021. SARS-CoV-2 cDNA was amplified using the tiling PCR approach developed by the ARTIC network [9] and sequenced on Oxford Nanopore sequencing devices (Figure 1C, methods section).

FIGURE 1 NEAR HERE

We first examined whether quantitative information on the viral load of the population can be deduced from wastewater sequencing data. Since SARS-CoV-2 enters wastewater essentially via human feces, co-amplification of a virus present in human feces appeared a promising strategy to obtain fecal load normalized SARS-CoV-2 data. We chose PMMoV (pepper mild mottle virus), a virus that is ingested with pepper-containing food. PMMoV is found at high concentrations in human feces and has previously been used for fecal load normalization of wastewater samples [START_REF] Bivins | Wastewater-Based Epidemiology: Global Collaborative to Maximize Contributions in the Fight Against COVID-19[END_REF][START_REF] Zhang | RNA viral community in human feces: prevalence of plant pathogenic viruses[END_REF]. To generate PMMoV normalized data, we included a PMMoV specific primer pair into pool 1 of the ARTIC SARS-CoV-2 primer panel and calculated for each sample the ratio between SARS-CoV-2 and PMMoV reads. Correlation of the SARS-CoV-2/PMMoV ratios with N-gene qPCR data for the individual samples revealed a batch variation between different sampling dates (Figure 2A, inset). We speculated that seasonal variations of pepper consumption or changes in food habits caused by COVID-19 related lockdowns might have led to fluctuations of fecal PMMoV load during the sampling period.

To correct for those batch variations we used qPCR data to calculate a correction factor for each sampling date (see Supplementary Table 2 for details). Figure 2A shows that there is a linear relationship between normalized SARS-CoV-2/PMMoV ratios and the number of viral copies derived from RT-qPCR. In this setting, PMMoV allows normalization of samples from the same date while seasonal variations require correction with qPCR data.

FIGURE 2 NEAR HERE

The weather conditions and the biological oxygen demand at 5 days (BOD5) of the wastewater remained stable at the different sampling times, as indicated by the stability of the flow rate and of the population equivalent (Supplementary Table 3). The observed variations of virus concentration in wastewater can therefore be interpreted as an appropriate indicator of the SARS-CoV-2 shedding rate by infected communities. Figure 2B indicates for the WWTP "Haliotis'' an increase of the virus concentration from October to February followed by a decline in March, which was consistent with other observations [12]. We noticed however different time courses and magnitudes of virus titers in different neighborhoods, possibly due to different modes of circulation of the virus. For instance, a strong peak appeared in "Les Moulins'' between December and February, and a single peak was observed in January and February in "Las Planas". These observations are consistent with the view that a higher population density, a characteristic of these socially deprived districts, is often associated with a higher circulation of the virus.

Two main objectives of our study were also : (1) to define the relative abundance of SARS-CoV-2 lineages and (2) to identify novel SARS-CoV-2 variants in wastewater samples.

Available bioinformatics workflows for SARS-CoV-2 sequencing data were not directly suitable for wastewater data since they were designed for the analysis of individual patient data, where each sample is assigned to one single lineage. We therefore designed a dedicated bioinformatics workflow where we first generated PCR primer trimmed, SARS-CoV-2 aligned reads with the ARTIC pipeline. We subsequently counted the mutation frequency at each position of the virus with the iVar software and assigned lineages based on mutational signatures defined in the Pangolin Github repository [START_REF] Rambaut | A dynamic nomenclature proposal for SARS-CoV-2 lineages to assist genomic epidemiology[END_REF] with a custom pipeline written in R (detailed in the methods section). The associated code is available at https://github.com/ucagenomix/cagablea. A global overview of the mutations found in the samples and the associated lineages is presented in Figure 3.

The percentage of a lineage was defined as the median of the ratios (number of reads)variant / (number of reads)total for the mutations defining the variant (Figure 4A). When mutations defining a lineage were located on amplicons derived from two distinct primer pools (distinct PCR amplifications, Figure 1), the lineage assignment correlated well when the data for each primer pool was analyzed separately. In some cases, a variant was defined by several mutations located on the same amplicon (B.1.351: Figure 4C; B.1.1.7: Supplementary Figure 1A). The search for multiple mutations on the same read allowed us to detect weakly expressed lineages with a high degree of confidence (B.1.351, Figure 4C). The available sequencing information was generally sufficient to ensure robust lineage determination (Figure 3, Figure 4B-C). Nine lineages were enough to cover most of the diversity found in Nice (Figure 4A). From October 2020 to the end of January 2021, the most abundant lineage sporadically in some neighborhoods, but in general, their frequency subsequently decreased rapidly (Figure 3). Most striking was the transient appearance of the A.27 lineage in the "Bon Voyage" neighborhood in January 2021, where it represented half of the virus population (Figure 3 and Figure 4). This lineage had been previously reported in the French overseas department Mayotte, then in Créteil [START_REF] Fourati | Novel SARS-CoV-2 Variant Derived from Clade 19B, France[END_REF] and the detection of this variant in Nice in January may be related to the arrival of travelers from these areas for Christmas. A.27 totally disappeared in February in this neighborhood when B.1.1.7 became prominent. The cumulative fraction of the different lineages shown in Figure 4A was usually around 100%.

The missing signal usually corresponded to the B.1 or B.1.1l lineages which were not represented in Figure 4A, because these 2 variants overlap with most of the other observed lineages (Figure 3).

FIGURES 3 & 4 NEAR HERE

The most striking observation is the appearance of B.1.1.7 in Nice in early 2021. This lineage was found in January 2021 in "Les Moulins" where it represented 50% of the reads (Figure 3, Figure 4A,B). B.1.1.7 was then also detected in "West Nice" and "Haliotis" where it represented 10-20% of the virus population. This is consistent with the fact that these two wastewater networks are collecting sewage from "Les Moulins" (Figure 1A). B.1.1.7 was also detected at lower levels in "Carras", "Nice Etoile", "Magnan" and "West Gambetta". The fraction of this strain then strongly increased in the whole city in February to 81±1% of the total sequences. B.1.1.7 remained predominant in March, representing 86±2% of the viruses circulating in the territory.

Interestingly, the majority of the B.1.1.7 variant that appeared in January in Nice had an additional A522S (G23126T) mutation in the spike protein (Figure 4D). This A522S spike protein mutation had been reported in only 226 out of 27,268 (0.8%) and 207 out of 385,840 (0.05%) GISAID submitted B.1.1.7 sequences for France and the UK, respectively [START_REF] Elbe | Data, disease and diplomacy: GISAID's innovative contribution to global health[END_REF]. This variant remained predominant in February and subsequently declined in March (Figure 4D).

The initial prevalence of this elsewhere rare B.1.1.7 mutant in Nice suggests that the B.1.1.7 variant was introduced by one or a few carriers and subsequently spread initially in "Les Moulins" and later in the entire city. This variant then became gradually diluted by B.1.1.7 lacking the A522S mutation which is predominant in France. The decline of the A522S mutation in March also suggests that this additional spike mutation does not provide improved fitness when compared to the original B.1.1.7 lineage. This is in line with a recent SARS-CoV-2 / ACE2 interaction modelling study that only found a slight increase of SARS-CoV-2 receptor affinity for the A522S mutation [START_REF] Wang | Vaccine-escape and fast-growing mutations in the United Kingdom, the United States, Singapore, Spain, India, and other COVID-19-devastated countries[END_REF].

By contrast, the circulation of other VOCs, such as B.1.351 and P.1, remained much lower (Figure 4B, Supplementary Figure 2A-B). Figure 3 shows that some mutations that are associated with these lineages could be detected at the end of 2020. As we only detected a fraction of the mutations defining those VOCs, we concluded that they were not specific to VOCs at this time. In only a few cases, a consistent signature with enough VOC-specific mutations could be identified. This was for instance the case for "B.1.351" in "Carras", "East Nice" and "Harbor" in February 2021. This presence was further demonstrated by the detection of two B.1.351 specific mutations on the same reads (Figure 4C). The P.1 lineage represented about 10% of the reads in "La Madeleine" in February 2021, but it was not detected there in March anymore (Fig. 4B).

We finally compared the lineage composition of the October wastewater sample with those of matched clinical samples. The 1,481 swab samples from patients living in Nice that were tested by Biogroup during a period corresponding to our October wastewater sampling (between Oct 19 th and Oct 23 rd , Week 43) were analyzed. There were 1,327 negative and 154 positive SARS-CoV-2 RT-qPCRs. We sequenced 81 PCR positive samples that could be unambiguously assigned to one of the wastewater catchment areas and determined their lineages (Figure 4E, Supplementary Table 4). We found similar percentages for the four most abundant lineages in wastewater and in clinical samples that were collected in the same week in the same area.

In summary, we defined the major lineages of SARS-CoV-2 that circulated from October 2020 to March 2021 in Nice. We were able to identify the switch that occurred between January and February 2021, characterized by a rapid onset of the B.1.1.7 lineage, and a joint decrease of the previously dominant lineages. As of February 2021, B.1.1.7 became the major lineage present in the city, representing more than 80% of all sequences. Our work shows small, localized spikes of the variants of concern B.1.351 and P.1 which did not spread widely into the whole city.

Altogether, our results demonstrate the reliability of wastewater RNA sequencing and the interest of this approach for the monitoring of the emergence and spread of SARS-CoV-2 variants. A limitation of our study is that our analyses were done with RNA from 1 milliliter of wastewater. An increase of the amount of RNA led to RT/PCR inhibition likely due to the co-purification of inhibitors present in wastewater during virus concentration and RNA preparation. Future SARS-CoV-2 concentration techniques and RT/PCR inhibitor removal from the RNA preparations will certainly improve the sensitivity of the measurement, as previously done for detecting poliovirus [START_REF] Berchenko | Estimation of polio infection prevalence from environmental surveillance data[END_REF]. This approach was however sufficient in the present study because the concentration of the virus was high at that time in Nice. The results were robust, and we believe that internal controls such as PMMoV will contribute to the development of more quantitative approaches.

In conclusion, our results strongly position RNA sequencing of wastewater samples as a credible tool for analyzing the spread and evolution of SARS-CoV-2. Viral monitoring in wastewater is clearly not limited to the surveillance of the present SARS-CoV-2 pandemic.

The surveillance approach that we have set up can easily be transferred to other epidemiological programs, targeting future emerging viruses of concern. As such, it appears as a valuable monitoring tool to assess the effectiveness of sanitary measures aiming to curb virus transmission. Last, our work illustrates how the situation can vary between different neighborhoods of the same city. The raise of B.1.1.7 in "Les Moulins" in January coincides with a raise of A.27 in "Bon Voyage", and of B.1.474 in "Las Planas". The emergence of those particular lineages in three socially deprived neighborhoods occurred immediately after Christmas, a time of very active traveling. This possibly highlights the link existing between socio-economic inequalities between different neighborhoods and the spread of the virus [START_REF] Horton | Offline: COVID-19 is not a pandemic[END_REF]. In this perspective, we believe that wastewater epidemiology could in the future represent a particularly useful instrument to better face similar situations. For each sampling date, a correcting ratio was defined by the weighted average signal between the different sampling points in the wastewater treatment plant (detailed in Supplementary Table 2). The inset shows the relationship between the SARS-CoV-2 / PMMoV ratio and the average RTqPCR signal for N gene before normalization. (B) N gene RT-qPCR Cq values were used to assess the concentrations of virus in the WWTP (Haliotis, Nice wastewater treatment plant) and the different neighborhoods. Values are provided in copies of genomes per ml of wastewater (cp/ml). Results are shown for 4 areas. Full results are provided in Supplementary Table 2. The population in each neighborhood is indicated under each name in italic.

Acknowledgements

  was B.1.160. Additional lineages (B.1.177, B.1.367, B.1.474, B.1.221, A.27) appeared

  Figure 2

Figure 2 :

 2 Figure 2: Comparison between RTqPCR and sequencing quantifications. (A) Relationship between normalized SARS-CoV-2 / PMMoV ratio and the average RTqPCR signal for the N gene.For each sampling date, a correcting ratio was defined by the weighted average signal between the different sampling points in the wastewater treatment plant (detailed in Supplementary Table2). The inset shows the relationship between the SARS-CoV-2 / PMMoV ratio and the average RTqPCR signal for N gene before normalization. (B) N gene RT-qPCR Cq values were used to assess the concentrations of virus in the WWTP (Haliotis, Nice wastewater treatment plant) and the different neighborhoods. Values are provided in copies of genomes per ml of wastewater (cp/ml). Results are shown for 4 areas. Full results are provided in Supplementary Table2. The population in each neighborhood is indicated under each name in italic.

Data and code availability

All relevant data are accessible through a dedicated webpage. They have also been deposited in Sequence Read Archive under BioProject ID:PRJNA750263 which is accessible at: http://www.ncbi.nlm.nih.gov/bioproject/750263. All custom software, R figures and analysis scripts are available on Github (https://github.com/ucagenomix/cagablea).

Materials and Methods

Wastewater sample collection. The characteristics of the different sampling sites are detailed in Figure 1 and Supplementary Table 1, which shows the sizes of the different populations, flow rates and volumes of water for each sampling. Fourteen sampling sites were equipped with automatic sampling devices (Sigma SD900, Hach Company, Loveland, Colorado). The instruments were calibrated to perform twelve samplings per hour for 24 hours. In some cases, especially at the beginning of this study, a manual collection was performed (Supp. Table1). All samples were processed within 24 hours. Sixty milliliters of 0.22µm-filtered wastewater were ultracentrifuged (1 hour, 158,000g, 4°c; Figure 1C). Pellets were resuspended, in 200µl of water, incubated for 5 minutes at room temperature and mixed carefully by pipetting up and down. Viral RNA was extracted using the AllPrep PowerViral DNA/RNA kit (Qiagen, Hilden, Germany, #28000-50) or the RNeasy PowerMicrobiome kit (Qiagen, #26000-50) following the manufacturer's instructions. RNA was eluted in 100µl of RNase free water and stored at -80°C. All wastewater and clinical samples were handled in a biosafety level 2 laboratory under a biological safety cabinet using adapted personal protective and respiratory equipment. RNA extraction. Viral RNA from 81 nasopharyngeal swabs were extracted using the QIAamp Viral RNA minikit (Qiagen, #52904) according to manufacturer instructions.

Samples RNA were eluted with 60µl of AVE buffer from the kit (0.04% of sodium azide in nuclease free water) and stored at -80°C before further analysis.

Gene N1 RTqPCR.

For SARS-CoV-2 RNA quantification by RT-qPCR, the N1 region of the nucleocapsid gene was targeted using the primers and probes as well as the PCR conditions described in the CDC diagnostic panel (2020). RT-qPCR quantifications were performed using the qScript XLT One-Step RT-qPCR ToughMix (Quantabio) on CFX96 touch Real-Time PCR Detection System. The N1 probe was labelled at the 5'-end with the FAM reporter and with the BHQ-1 quencher at the 3'-end. For each sample, RT-qPCR analysis was systematically performed in duplicate of which the mean is reported in the study.

Prior to the RT-qPCR, extracted RNAs were diluted with RNase free water to bring the initial RNA amount back in a range of [50-150 ng] RNA quantity per reaction. When the RNA concentration was below 50 ng, the sample was by default diluted by a factor of 2, in order to dilute matrix-related inhibitors. An inhibition control of the RTqPCR reaction was systematically performed for each extracted sample. This control consisted of a dosed addition (~4000 copies/reaction) of the positive control EDX SARS-CoV-2 (Exact Diagnostic, supplied by Bio Rad) for each dilution analyzed. The amplification was considered as free of inhibition with the Cq 24 with a tolerance of 2 Cq. A standard curve, from 10 6 copies/μL to 1 copy/μL, was generated for every analysed plate using a control plasmid, N-nCoV-control-Plasmid (Eurofins). Standard curves were performed by serial 10 fold dilution with 1X Tris -EDTA. A negative control was analyzed on each plate, using SARS-CoV-2 Negative Run Control (Exact Diagnostic, supplied by Bio Rad). This control is formulated in a synthetic matrix and contains 75,000 copies/mL of human genomic DNA.

SARS-COV-2 genome amplification and sequencing.

Library preparation was performed using the nCoV-2019 sequencing protocol v3 LoCost [19] with a few modifications. For reverse-transcription, 8 µl of RNA and 2µl of LunaScript (NEB) was mixed and incubated at 25°C for 5 minutes followed by an incubation at 55°C for 10 minutes before holding at 4°C. The tables generated by iVar were filtered for positions that were covered by at least 100 reads and for alterations with an adj p-value (Fisher's test) below 10 -8 .

A custom R-script was implemented to fix a known issue with the iVar output, the incorrect assignment of deletion frequencies and to code the mutations events in amino acid format.

A database with the mutations associated with known lineages was built from the publicly available data [21].

Every mutational event in the sample was associated with mutations described in the database and matched to one or several lineages. A list of specific mutations for each lineage of interest is provided in Supplementary Figure 3. For each of these positions we calculated the ratio of the number of mutant reads to the total number of reads. The frequency of each lineage was then defined as the median of these lineage-related ratios. We usually noticed that the total percentage of the different lineages was close to 100%. Percentages below 100% were probably due to the presence of additional lineages that could not be properly identified ( e.g. lineages that lack characteristic mutations that are not shared with other lineages). We also noticed in some cases total percentages that were slightly above 100%. This is probably due to experimental noise of the measurements. As this effect was moderate, we did not correct this signal, to avoid overcorrection and to keep data processing and interpretation simple. October November December January February March

Supplementary Figures 3BC indicate that

Q C P F T N E S M L H D G I I A A Q C P F T N E S M L H D G I A K J I Q C P T F N S E H D B G M A L K U J I Q C P O O T N E S M L H D G U I I A K J B I A B J Q R C O P F T N E S L M H D U G K I A B J Q R C O P F T N E S L M H D K U G I

A 4.