Mohammed Nadhir

D Cherief

Fahmi Zairi
email: fahmi.zairi@polytech-lille.fr

Ning Ding

Jean-Michel Gloaguen

Moussa Naït-Abdelaziz

Mohamed Benguediab

Fahmi Zaïri

Plasticity and thermally-induced recovery in polycarbonate

Keywords: Amorphous thermoplastics, Thermally-induced recovery, Viscoplasticity, Viscohyperelasticity

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

Plasticity is a common term to designate the propensity of a material to undergo unrecoverable alteration of the microstructure by the action of a mechanical loading: the material remains deformed after the mechanical loading is removed. Although amorphous polymers exhibit what appears as a plastic deformation well described now in the literature, paradoxically, the "apparent" plasticity may be regarded as a frozen process only at the loading time-scale since it may be fully or partially recovered in response to a suitable stimulus such a temperature change. The purpose of this contribution is to come back to this 2 subject by studying the thermally-induced recovery of the apparent plasticity, in the specific case of polycarbonate, by an approach combining constitutive modeling and experiments. This strain recovery is now commonly termed "shape-memory" phenomenon and was firstly introduced in the mid of the last century. It is only from the past twenty years that the research on shape-memory polymers has seen an explosive growth in both academics and industrial fields [START_REF] Lendlein | Shape-memory polymers[END_REF][START_REF] Hu | Recent advances in shape-memory polymers: Structure, mechanism, functionality, modeling and applications[END_REF][START_REF] Zhao | Recent progress in shape memory polymer: New behavior, enabling materials, and mechanistic understanding[END_REF]. They are merged as a new class of self-repairing and intelligent materials with a very wide range of potential application areas from biomedical engineering for minimal invasive surgery to aerospace industry. Nonetheless, all polymers have a thermally-induced shape-memory capacity even under large-strains. The shape-memory behavior of a polymer is due to the flexibility of the polymer chains and their capacity to store and release mechanisms through treatments around a defined transition temperature. For amorphous networks, the transition temperature is the glass transition temperature Tg. The thermally-activated shape-memory effect may be evidenced using conventional techniques as visualized in Figure 1 when the polymer specimen passes through four steps to form a cycle with different temporary states from the initial state I to the recovered state V1 . The compression is chosen because of the homogeneous nature of this deformation mode within large deformation behavior. Below Tg, the glassy specimen is compressed at a given strain-rate ε& to the desired strain (state II) in which the post-yield response is strongly nonlinear and generally governed by two features: a strain-softening immediately after yield-peak followed by a strain-hardening at large-strains. The specimen is then unloaded down to zero-stress 0 σ = at the same absolute strain-rate ε& (state III). Upon unloading, a part of the stored energy is released. The nonlinearity in unloading is indicative of the instantaneous recovery of an amount of the inelastic strain beyond that associated with the linear elasticity. At the end of the unloading process, a substantial remanent strain is retained which necessitates a very long recovery time as long as the temperature is held constant. To activate the recovery process, the deformed specimen is firstly cooled down to room temperature T0 (state IV) in a mechanical loading free condition and then reheated above Tg under zero-stress (state V). The recovery process is strongly dependent on the previous thermo-mechanical history (temperature, strain-level and strainrate) and the heating conditions in terms of recovery temperature and heating rate. The constitutive model of [START_REF] Haward | The use of a mathematical model to describe isothermal stress-strain curves in glassy thermoplastics[END_REF], extended later by [START_REF] Boyce | Large inelastic deformation of glassy polymers. Part I: Rate dependent constitutive model[END_REF] into a 3-D finite strain formulation, appears to be the first effort to formulate a constitutive theory for amorphous thermoplastics below their glass transition temperature and has been early developed based on the observation of the recovery. It has been the inspiration of a multitude of constitutive models for the time-temperature dependence of the yielding behavior in continuum terms invoking mechanical analogs by using combinations of elements, usually linear elastic springs, viscoplastic dashpots and nonlinear elastic springs [START_REF] Arruda | Effects of strain rate, temperature and thermomechanical coupling on the finite strain deformation of glassy polymers[END_REF][START_REF] Buckley | Glass-rubber constitutive model for amorphous polymers near the glass transition[END_REF][START_REF] Tervoort | A constitutive equation for the elasto-viscoplastic deformation of glassy polymers[END_REF][START_REF] Boyce | Constitutive model for the finite deformation stress-strain behavior of poly(ethylene terephthalate) above the glass transition[END_REF][START_REF] Lion | Constitutive modelling in finite thermoviscoplasticity: A physical approach based on nonlinear rheological models[END_REF][START_REF] Richeton | Modeling and validation of the large deformation inelastic response of amorphous polymers over a wide range of temperatures and strain rates[END_REF][START_REF] Ames | A thermo-mechanically coupled theory for large deformations of amorphous polymers. Part II: Applications[END_REF][START_REF] Anand | A thermo-mechanically coupled for large deformations of amorphous polymers. Part I: Formulation[END_REF][START_REF] Belbachir | Modelling of photodegradation effect on elasticviscoplastic behaviour of amorphous polylactic acid films[END_REF]Srivastava et al., 2010a;[START_REF] Zaïri | Constitutive modelling of the large inelastic deformation behaviour of rubber-toughened poly(methyl methacrylate): Effects of strain rate, temperature and rubber-phase volume fraction[END_REF][START_REF] Zaïri | A physically-based constitutive model for anisotropic damage in rubber-toughened glassy polymers during finite deformation[END_REF][START_REF] Bouvard | An internal state variable material model for predicting the time, thermomechanical, and stress state dependence of amorphous glassy polymers under large deformation[END_REF][START_REF] Bouaksa | Molecular chain orientation in polycarbonate during equal channel angular extrusion: Experiments and simulations[END_REF][START_REF] Nguyen | A large strain hyperelastic viscoelastic-viscoplastic-damage constitutive model based on a multimechanism non-local damage continuum for amorphous glassy polymers[END_REF][START_REF] Gudimetla | A finite strain thermodynamically-based constitutive framework coupling viscoelasticity and viscoplasticity with application to glassy polymers[END_REF][START_REF] Mahjoubi | A micro-macro constitutive model for straininducedmolecularordering in biopolymers: Application to polylactide over a wide range of temperatures[END_REF][START_REF] Mahjoubi | Strain-induced phase transformation in poly(lacticacid) across the glass transition: Constitutive model and identification[END_REF]. As a common point, a hyperelastic strain-hardening resistance is associated to an elastoviscoplastic formulation to account for the intermolecular resistance. Mathematically, the connection between both resistances vary from a model to another and is realized either by considering the orientation hardening as a kinematic back-stress or as an additive stress. Physically, the models consider generally that the orientation hardening originates from the rubber-elastic response of the macromolecular chains between junctions in the macromolecular network that persist through a temperature range, called physical entanglements in the case of an amorphous thermoplastic and frozen in when it is cooled from the melt state. The concept of physical entanglements and rubber-elastic analogy in glassy polymers is now commonly accepted to explain the strain-hardening at large-strains. The strongest argument is the increase of the strainhardening with the entangled polymer network density due to enhanced straightening out of chains between entanglements and associated orientation of the macromolecular network. The stress contribution of the orientation hardening is generally constitutively described using the rubber-elastic network theories where the cross-links play the role of junctions in the macromolecular network. The effort of modeling for amorphous thermoplastics continues to date, but it is generally focused on the reproduction of the stress-strain behavior, sometimes spanning the glass transition temperature, but generally without being interested in the thermally-induced shape-memory capacity. The earlier modeling attempts to describe the thermally-activated recovery effects are due to [START_REF] Tobushi | Thermomechanical constitutive modeling in shape memory polymer of polyurethane series[END_REF][START_REF] Tobushi | Thermomechanical constitutive model of shape memory polymer[END_REF] and [START_REF] Lin | Shape-memorized crosslinked ester-type polyurethane and its mechanical viscoelastic model[END_REF] studying shape-memory polyurethane based-polymers in the framework of thermo-viscoelasticity combined with a frictional element. Later, two main classes of models have been proposed to capture the thermallyinduced shape-memory behavior of amorphous polymers considering either the viscosity temperature-change approach [START_REF] Diani | Finite strain 3D thermoviscoelastic constitutive model for shape memory polymers[END_REF][START_REF] Nguyen | A thermoviscoelastic model for amorphous shape memory polymers: Incorporating structural and stress relaxation[END_REF][START_REF] Castro | Effects of thermal rates on the thermomechanical behaviors of amorphous shape memory polymers[END_REF]Srivastava et al., 2010b;[START_REF] Chen | Influence of thermoviscoelastic properties and loading conditions on the recovery performance of shape memory polymers[END_REF][START_REF] Su | A 3D finite strain viscoelastic constitutive model for thermally induced shape memory polymers based on energy decomposition[END_REF][START_REF] Mao | A viscoelastic model for hydrothermally activated malleable covalent network polymer and its application in shape memory analysis[END_REF] or the rubbery-glassy phase-transition approach [START_REF] Liu | Thermomechanics of shape memory polymers: Uniaxial experiments and constitutive modeling[END_REF]Chen andLagoudas, 2008a, 2008b;[START_REF] Qi | Finite deformation thermomechanical behavior of thermally-induced shape memory polymers[END_REF]Ghosh andSrinivasa, 2013, 2014;[START_REF] Boatti | A three-dimensional finite-strain phenomenological model for shape-memory polymers: Formulation, numerical simulations, and comparison with experimental data[END_REF][START_REF] Park | Three-dimensional constitutive model for shape memory polymers using multiplicative decomposition of the deformation gradient and shape memory strains[END_REF][START_REF] Yang | Temperature and rate dependent thermomechanical modeling of shape memory polymers with physics based phase evolution law[END_REF][START_REF] Li | A viscoelastic constitutive model for shape memory polymers based on multiplicative decompositions of the deformation gradient[END_REF]. The first class of models introduces the temperaturedependence of the molecular mobility and the relaxation time. The second class of models is the widely used approach due to its simplicity. It introduces a phase-transition hypothesis considering the amorphous network as a continuum mixture of a glassy phase and a rubbery phase using the concept of phase volume fraction to describe the evolution of the microstructure. The glassy and rubbery phases are stable below and above the transition temperature, respectively, and the transitions between the two superimposed phases with a temperature change provide a description of the strain storage and recovery during a shapememory cycle. The phase-transition approach, originally developed for shape memory alloys, finds a strong physical meaning when the shape-memory behavior concerns polymer medium with a phase transformation, i.e. crystallizable polymers (Guo et al., 2018a;[START_REF] Moon | A thermodynamic framework for the modeling of crystallizable triple shape memory polymers[END_REF]. In an original amorphous medium that remains amorphous, the phase-transition approach remains conceptual and the lack of physical meaning constitutes the main limitation. Most of models deal with the thermally-induced shape-memory of thermoplastic elastomers or thermoset polymers and are limited to certain classes of constitutive responses. Our approach differs from the above models. With the aim of developing an approach that is as general as possible and physically-based, the presented elastoviscoplastic-viscohyperelastic model includes the variation of the intermolecular and network properties with temperature and their abrupt change when the temperature traverses the glass transition temperature. The constitutive equations combine existing elements from several theories (elasticity, viscoplasticity and viscohyperelasticity) that are presented in a unifying manner and in a new way according to the model structure in order to facilitate finding the key to the thermallyactivated strain recovery in connection with the microstructure and to give a clearer definition of the "plasticity" in amorphous polymer systems. Our approach is simplest than previous models and does not need to consider sub-networks and the introduction of an additional mathematical formulation for strain recovery during heating. While our approach remains enough simple, it is based upon a physical understanding of the material behavior by considering the underlying thermo-mechanical mechanisms and their abrupt change in the temperature dependence when the material changes from the glassy state to the rubbery state. For the development of the thermo-mechanical model, compression experiments on a polycarbonate were conducted over a wide range of temperatures across the transition temperature. The model prediction is studied for zero-stress creep recovery. The present paper is organized as follows. In Section 2, the finite strain constitutive theory is presented in the nonlinear continuum thermo-mechanical framework2 . In Section 3, the attention is devoted to the quality of the model fit to the experimental stress-strain curves and to the predictive model capability for recovery response. We close in Section 4 with final remarks.

Theory

The apparent plasticity is regarded as a thermally-activated process determined by elementary molecular events related to intermolecular and intramolecular motions, the thermallyactivated recovery being made possible by the introduction of their variation with temperature. Figure 2 provides a schematic view of an amorphous network along with the elementary molecular events occurring inside. The restriction imposed on molecular chain motion due to neighboring chains leads to intermolecular interactions resulting in a first resistance to small-moderate deformations. The latter resistance (intermolecular resistance I) is responsible for the initial stiffness and the yield events. At larger deformations, the orientation and relaxation processes of the macromolecular network lead to a second resistance (molecular network resistance N) responsible for the remarkable strain-hardening response.

Molecular events and properties

It is now well agreed that the glass transition is the key physical feature for the shape-memory effects of amorphous networks. Nonetheless, this subsection introduces fundamental relations between all elementary molecular events and properties for the sake of completeness in the description of polymer deformation. The nature of the molecular mobility in amorphous thermoplastics results from intermolecular barriers to chain-segment rotation of two distinct molecular sub-segments related to thermally-activated α and β processes acting at two different levels3 . The activation energy required for the very local movement of the β process is low, and its relaxation is activated before the α process. The two relaxations are key features for the polymer deformation and a physically consistent description is needed. An Arrhenius-type expression describes the strain-rate dependence of the β transition temperature (Richeton et al., 2005b):

( ) The α transition is also referred to as the glass transition which follows the Williams-Landel-Ferry (WLF) free-volume theory (Richeton et al., 2005b):

1 1 ln ref ref k T T H β β β ε ε ε   = +   ∆   & & & (
( )

1 2 1 log log ref ref g gref g g T T c c ε ε ε ε ε -       = - +           
in which gref T is the α transition temperature at ref ε& and, 1g c and 2 g c are the WLF coefficients.

To complete the description of the temperature variation, the melting temperature m T is described by a phenomenological formula (Richeton et al., 2005b):

( ) 2 1 10 log m mref ref T T ε ε ε -     = +             & & & (3)
in which mref T is the melting transition temperature at ref ε& .

The elasticity and the plastic yielding of amorphous thermoplastics are strongly affected by strain-rate and temperature, with a similar dependency. As shown in Figures 3a and3b, these two quantities decrease with temperature and increase with strain-rate. The latter effects must be correlated to the transition temperatures. The variation of the Young's modulus E with temperature θ and strain-rate ε& is given by the following expression, valid from fully glassy state to fully rubbery state (Richeton et al., 2005b):

( ) ( ) ( ) ( ) , , , , g m T T T E E E E β θ ε θ ε θ ε θ ε
the following formulae using Weibull statistics (Richeton et al., 2005b):

( ) ( ) ( ) ( ) ( ) 1 1 2 , exp β β θ θ ε ε ε ε       = - -          & & & & m T E E E T (5) ( ) ( ) ( ) ( ) ( ) 2 2 3 , exp θ θ ε ε ε ε       = - -          θ θ ε ε ε       = -          & & & m m T m E E T (7)
where 1 m , 2 m and 3 m are the Weibull shape parameters of each transition region and, ( )

1 E ε& , ( ) 2 E ε& and ( ) 3
E ε& are the instantaneous moduli at the beginning of each transition region (Richeton et al., 2005b):

( )

1 1 1 log ε ε ε     = +             & & & ref ref E E s (8) ( ) 2 2 1 log ε ε ε     = +             & & & ref ref E E s (9) ( ) 3 3 1 log ε ε ε     = +             & & & ref ref E E s (10)
where 1ref E , 2ref E and 3ref E are the instantaneous moduli at a chosen reference strain-rate ref ε& and, s is a parameter characterizing the strain-rate sensitivity of the stiffness modulus assumed to be the same for the three transition regions. Figure 3a presents the Young's modulus variation with temperature of polycarbonate, using the above equations, from fully glassy state to fully rubbery state. In the isotropic case, the elasticity may be completely defined by the Young's modulus E and the Poisson's ratio ν .

The latter material property was assumed to be strain-rate independent but temperature dependent4 :

( ) ( ) ( ) 2 2 10.3 0.36 0.139exp 10.3 θ θ ν θ   - + -   = +   +     gref gref T T for g T θ < (11)
below g T , and

( ) 0.499 ν θ = for g T θ ≥ (12)
above g T . Once the energy barrier to molecular mobility is overcome the elastic behavior gives way to the plastic yielding through the thermally-activated local rotation of molecular segments. Figure 3b presents the yield stress variation with temperature of polycarbonate from fully glassy state to fully rubbery state. At low temperatures, the increase of the yield strength with strain-rate is correlated to the decrease of the molecular mobility leading to the stiffening of the chains. This increase of the yield strength is associated to the activation of both α and β processes. However, at high temperatures and low strain-rates, the intermolecular resistance is mainly due to the α process while the β process is hind. The physics of plastic yielding is described by the accumulated viscoplastic strain rate vp I γ& that differs on both sides of g T :

(

)

vp vp int ˆ, , I I H V β γ γ σ = ∆ & & for g T θ < (13)
below g T , and

( )

vp vp ˆ, I I H V β γ γ = ∆ & & for g T θ ≥ (14) above g T .
The term H β ∆ is again the β activation energy, V is an activation volume and int σ is the internal intermolecular stress frozen in the deformed network in the glassy state (see Figure 2). The internal stress evolves from its initial value

( ) int =0 t σ
to the value ps σ referring to the preferred structural state of the polymer chains according to the following evolution law:

vp int int 1 I ps h σ σ γ σ   = -       & & (15)
The evolution ( 15) is introduced to account for the strain-softening response in the glassy state manifested by a drop in stress with strain upon initiation of the plastic flow. The term h 7 is the strain-softening slope. The ratio

( ) int 0 ps t σ σ =
is assumed constant and the temperature dependence of ( ) int 0 t σ = is given by a linear law [START_REF] Richeton | Influence of temperature and strain rate on the mechanical behavior of three amorphous polymers: Characterization and modeling of the compressive yield stress[END_REF]:

( ) ( ) int int 0 0 Im t m θ σ σ θ ασ = = = -+ (16)
where

( )
int 0 θ σ = is the internal intermolecular stress at the temperature θ = 0 K, m is the sensitivity parameter vis-à-vis the temperature and α is the sensitivity parameter vis-à-vis the intermolecular mean stress Im σ .

At larger deformations, the molecular network experiences an orientation process leading to the rapid stress increase at the macroscopic scale. The latter material feature is directly related to the limiting extensibility of the molecular network and can be identified by two physical ingredients: the initial hardening modulus r c

C n kθ = 17) and ( 18) for the orientation hardening but by introducing them in a kinematic back-stress. In addition, the entangled molecular network orientation process translated by the kinetics ( 17) and ( 18) enters in competition with the molecular network relaxation process due to chain reptation allowing decreasing the stress5 . The physics of molecular network relaxation is described by the accumulated viscous strain rate v N γ& :

( )

v v N N C γ γ = & & (19)
As presented in Figure 3c, the term C is a viscosity parameter that linearly increases with temperature with the same slope on both sides of the glass transition [START_REF] Zaïri | Constitutive modelling of the large inelastic deformation behaviour of rubber-toughened poly(methyl methacrylate): Effects of strain rate, temperature and rubber-phase volume fraction[END_REF]:

0 C C cθ = + (20)
where 0

C and c are constants.

Constitutive equations

This subsection presents constitutive relationships between the macro-stress response and the set of physically-interpretable material parameters introduced above along with their remarkable change when the temperature traverses the glass transition. The model is developed from the concept that the intermolecular resistance I and the molecular network resistance N participate to the total resistance to deformation in polymers by acting in parallel.

The total Cauchy stress σ is thus the sum of both contributions:

= + σ σ σ I N

(21) The two basic physical mechanisms are, respectively, formulated using the elastoviscoplastic theory, to reproduce the initial stiffness and the yield events, and the viscohyperelastic theory, to reproduce the remarkable strain-hardening at large-strains6 . The general form of the model is illustrated schematically in Figure 2 for a one-dimensional visualization analog of the polymer deformation model: the resistance I is constituted by a linear spring in series with a viscoplastic damper and, the resistance N is constituted by a hyperelastic spring in series with a viscous damper. The two resistances, acting in parallel, are joined in series to a thermal branch responsible for the stress-free thermal expansion. The corresponding kinematics of large-strain thermo-mechanics is recalled in Appendix A. ( )

( ) ( ) ( ) ( ) ( ) ( ) , 2 1 2 2 1 e I ik jl il jk ij kl ijkl E θ ε ν θ δ δ δ δ δ δ ν θ ν θ   = + +   - +     C & (23)
in which the term δ denotes the Kronecker-delta symbol.

The network Cauchy stress σ N is expressed as a function of the corresponding strain with hyperelastic spring analog. It is constitutively expressed using a rubber-elastic network formulation introducing the initial hardening modulus r C and the average number of rigid links rl N . The Arruda-Boyce 8-chain model is used [START_REF] Arruda | A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF]:

( ) ( ) ( ) ( ) ( ) h 2 1 h h h h 1 3 T rl N N N N N N N r r l C N J N θ λ λ θ θ λ -       = -         σ F F I L (24) in which ( ) h det = F N N J represents the network volume change, ( ) ( ) h h h tr 3 T N N N λ = F F
is the average hyperelastic network stretch on a chain where ( )

1 3 h h N N N J - = F
F and I is the identity tensor. The inverse Langevin function is given by a Padé approximation:

( ) ( ) ( ) For the accumulated viscoplastic strain rate vp I γ& , the mathematical formulation of [START_REF] Richeton | Influence of temperature and strain rate on the mechanical behavior of three amorphous polymers: Characterization and modeling of the compressive yield stress[END_REF]Richeton et al. ( , 2005a) is retained:

1 2 2 3 1 x x x x - ≈ - - L .
( ) int vp 0 exp sinh 2 I n I V H R k β σ γ γ θ θ   - ∆   = -           σ & & for g T θ < (27)
below g T , and

( )

1 vp 0 2 ln10 exp exp sinh 2 g g n I I g g g c T H V RT c T k β θ γ γ θ θ     × - ∆     = -         + -       σ & & for g T θ ≥ (28)
above g T .

The term 0 γ& is the pre-exponential factor, R is the perfect gas constant and n is a parameter characterizing the cooperative movement of the chain segments. The time-temperature dependence of the orientation hardening at large-strains is constitutively described using the expression (24) for the entangled polymer network orientation process combined with a viscous formulation for the relaxation process. The accumulated viscous strain rate v N γ& is formulated within a viscoplastic theory with no yield surface [START_REF] Pyrz | Identification of viscoplastic parameters of phenomenological constitutive equations for polymers by deterministic and evolutionary approach[END_REF], the viscous deformation being non-zero from the beginning of the mechanical loading. Motivated by the concept of the reptational Brownian motion, the Bergstrom-Boyce power law expression for stress/strain-assisted molecular network relaxation rate is retained [START_REF] Bergstrom | Constitutive modeling of the large strain time-dependent behavior of elastomers[END_REF]:

( )

1 v v 1 N N N C γ θ λ - = - σ & (29)
where

( ) ( ) v v v tr 3 T N N N λ = F F
is the viscous network stretch with an initial value that deviates slightly from 1.0 in the computations to avoid singularity.

Thermally-activated strain recovery

Figure 4 presents the loading-unloading-recovery process applied to the specimen. At the glassy state (θ <Tg), the macromolecular chains have low mobility, so that the intermolecular resistance I is enhanced and the material stores a substantial remanent strain when the mechanical loading is removed as illustrated in Figure 1. The amount of this apparent plastic deformation depends on the temperature and the strain-rate. In this unloaded state, the elastic intermolecular deformation is totally recovered, i.e. e I = F I , and the internal intermolecular stress int σ is frozen in the material. An amount of the network deformation N F is also recovered (due to the parallelism of intermolecular and network resistances) which leads to the stress nonlinearity upon unloading, the remaining network deformation contributing to the remanent strain. When the polymer is cooled down to room temperature (θ =T0), it becomes in complete frozen state and the total stress σ becomes null. Then, Eq. ( 21) writes: 25) and ( 27):

I N + = σ σ 0 at 0 T θ = (
-1 In Eq. ( 30), the network stress σ N is given by Eq. ( 24). Inherent to its formulation, the model considers that the two resistances act in parallel and without interaction. Nonetheless, a correlation exists between intermolecular and network properties. Indeed, as shown in Figures 3a, 3b and 3c, the kinetics with temperature increase is due to increase in molecular mobility, chain orientation and relaxation process. At the rubbery state (θ >Tg), the macromolecular chains can undergo large random conformational motions which increase their mobility, such that the intermolecular resistance I becomes negligible. More precisely, the response is governed by the orientation/relaxation network resistance N. In this region, the macromolecular chains are dominated by entropy as indicated by the increase of r C in the right side of the glass transition. Also, the chain viscosity translated by C is important: in this region, the macromolecular chains slide more easily against each other and take the capability to significantly change conformation by Brownian motion in a combination of reptation motion, as considered in the Doi-Edwards concept and constitutively represented by Eq. ( 29).

1 vp vp vp int vp 0 2 2 arsinh exp n I I I I I H k V R β γ θ σ γ γ θ    ∆      = +             σ F F & & & & (
The shape-memory behavior is therefore a consequence of a remarkable change in thermomechanical properties at the chain-scale in response to a temperature change. In particular, during the heating process from glassy state up to rubbery state, the inelastic features change quickly in response to a temperature variation.

As the true stress-strain curves, the recovery response may be considered as an intrinsic feature obtained by the following zero-stress creep condition above glass transition:

I N + = σ σ 0 and g T θ > (32)
where the intermolecular stress I σ is again expressed by Eq. ( 31) and the network stress σ N by Eq. ( 24). The accumulated viscoplastic strain rate vp I γ& in Eq. ( 31) is governed by Eq. ( 28) in which I σ is replaced by N -σ according to the condition (32). Note that I N = σ σ . Via the material kinetics of Section 2.1 introduced in Eq. (32), the model makes possible to consider that the deformed macromolecular structure keeps "memory" of its original position during heating above glass transition, the network resistance being used as the driving stress that continuously activates the recovery process during zero-stress creep above glass transition. By this way, the polymer regains the entropy that was lost during the chain orientation and relaxes to an equilibrium configuration in which it recovers in a great extent the initial conformation by "memory" effect. To further illustrate our approach, the viscoplastic intermolecular stretch is plotted as a function of the viscous network stretch in Figure 4. The increase in plastic strain during isothermal mechanical loading ( I II → ) is followed by a frozen process during unloading ( II III → ). The heating above glass transition combined with zero-stress creep ( IV V → ) metamorphoses the apparent plasticity into a recoverable process while the viscous network stretch exhibits continuous creep behavior as a manifestation of the reptational Brownian motion.

Application to polycarbonate 3.1. Experiments

A series of experimental tests have been performed on polycarbonate cylinders to get the stress-strain and recovery response assessing the thermally-induced shape-memory cycle as already illustrated in Figure 1. The experiments were carried out on an electromechanical model-5567 Instron universal testing machine equipped with an oven. Cylindrical specimens of 8.00 mm length and 6.3 mm diameter were machined from polycarbonate cylinders purchased from Goodfellow7 . The length to-diameter ratio of the specimens allows getting the stress-strain response at large-strains by avoiding both buckling and non-homogeneous deformation. Lubricant was used between the compression platens and the specimen surfaces (smooth and parallel) to avoid both barreling and shearing phenomena during the compression. Three kinds of tests were performed:

• Loading-unloading tests at different strain-rates and temperatures: The specimens were uniaxially loaded at a constant strain-rate until a prescribed strain and then unloaded at the same absolute strain-rate until zero-stress. The tests were conducted at three strain-rates ε& (0.01 s -1 , 0.001 s -1 and 0.0001 s -1 ) over a wide range of compression temperatures θ : 298 K, 333 K, 363 K, 393 K and 433 K. This experiment reproduces the behavior of the states I, II and III in Figure 1. • Zero-stress creep recovery above glass transition: The specimens, previously compressed below g T , were cooled down to room temperature by natural cooling and then submitted

to creep experiments at zero-stress for two heating temperatures: 433 K ( g T + 10 K) and 473 K ( g T + 50 K). This experiment reproduces the behavior of the states IV and V in Figure 1. The recovery R will be defined as: max 1-100

ε ε   = ×     R (33)
where ε and max ε are, respectively, the remanent strain stored in the material and its maximum value.

• Successive constant-temperature loading-unloading tests below glass transition: The specimens were successively compressed until a prescribed strain at four temperatures: 298 K, 333 K, 363 K and 393 K. Between each loading-unloading test, the compressed specimen was cooled down to room temperature by natural cooling. The strain was small enough (10%) to avoid orientation hardening and the strain-rate was fixed to 0.001 s -1 . Two cases were considered: increasing compression temperatures from 298 K to 393 K and decreasing compression temperatures from 393 K to 298 K. In what follows, the temperature and strain-rate dependent loading-unloading results are used to obtain the kinetics of the model constants while the successive constant-temperature loading-unloading results and the thermally-activated recovery results serve to verify the model predictive capabilities.

Model fit

Figure 5 presents the deformation behavior of polycarbonate for different temperatures and strain-rates, up to a true strain-level of 70%. The polymer exhibits the expected dependencies on temperature with two distinct responses below and above the glass transition. A zoom is provided in Figure 5d for the highest compression temperature. Below the glass transition, the material exhibits the stress-strain behavior characteristics of a typical amorphous thermoplastic in the glassy state from small to large-strains. The strain-rate affects the plastic flow stress (yield-peak and strain-softening) but both stiffness modulus and hardening modulus present a relatively weak dependence on strain-rate. The plastic flow stress drops significantly upon a temperature increase. The remanent strain, for which stress is zero at unloading, is slightly affected by the temperature. When the temperature increases further and exceeds the glass transition, the polymer exhibits a very low yield strength and responds as a hysteretic-rubber-like material for which the material viscosity is still significant as evidenced by the strain-rate dependency. The remanent strain is significantly different as the strain-rate changes. Furthermore, it can be noticed that the remanent strain seems higher than that generally observed for cross-linked network systems such as rubbers (e.g. Guo et al., 2018d). The intermolecular and network constants are identified separately thanks to the model structure. The intermolecular properties are mainly defined by the stiffness modulus and the yield strength. They are plotted in Figures 3a and3b as a function of temperature in the form of symbols. The first quantity is deduced from the initial slope of the stress-strain curve. The second quantity is deduced from the zero slope of yield-peak below the glass transition and the linear deviation above the glass transition. The stiffness modulus remains almost constant in the glassy region whereas the yield stress decreases monotonically with temperature. Both quantities change abruptly when the temperature traverses the glass transition. Table 1 provides the values of the model constants. The glassy-rubbery stiffness variation with temperature allowed identifying the model constants in Eqs. ( 1)-( 10). A reference strain-rate ref ε& of 1 s -1 was used. The model constants in Eqs. ( 16), ( 27) and ( 28) were extracted by means of a least-square adjustment of the yield stress as a function of temperature and over a range of strain-rates. The model constants in Eq. ( 15) were determined by fitting the model to the strain-softening region of the stress-strain data in the glassy state. Although, they were assumed constant vis-à-vis temperature, the model accounts for the temperature dependence of the strain-softening phenomenon. This feature may be illustrated in Figure 3d in which the work hardening rate, defined as the slope of the stress-plastic strain curve, is plotted as a function of stress. Figure 3c presents the three network parameters as a function of temperature. The model constants in Eqs. ( 17), ( 18) and ( 20) were successively iterated in order to achieve optimal agreement between the model and the orientation hardening region on both sides of the glass transition. The thermal conductivity effect could have been taken into account as in (Guo et al., 2018b(Guo et al., , 2018c) ) but only homogeneous simulations under isothermal conditions were considered. Figure 5 shows the ability of the model to reproduce the main features of the stress-strain curves up to large-strains for the different compression temperatures and strain-rates. That indicates that the proposed approach provides a good approximation of the intrinsic elastoviscoplastic-viscohyperelastic response of this polymer and the basic underlying thermo-mechanical mechanisms. Quite interestingly, the model accounts for the abrupt change in intermolecular and network properties when the temperature traverses the glass transition temperature and the thermoplastic/rubber-like response transition. This point is fundamental in the aim reproducing the thermally-activated shape-memory effect.

Model prediction

The loading-unloading response, followed by cooling to T0, freezes the macromolecular structure and the specimen shape by the appearance of an apparent plastic deformation comprising the viscoplastic intermolecular deformation and the network deformation not released after unloading. This frozen state is a non-equilibrium configuration for which the molecular network relaxation towards equilibrium is extremely long: an extremely long recovery time is needed as long as the temperature is held constant. Successive compressions are presented in Figure 6 for increasing compression temperatures (from 298 K to 393 K) and decreasing compression temperatures (from 393 K to 298 K). Each constant-temperature curve is followed by a cooling to T0. The experiments were replicated by the model using the constants calibrated in Section 3.2 by considering two kinds of simulations under isothermal conditions: (i) in the first one, named intermittent simulations, the material is considered as a stress-free/strain-free medium at the beginning of each cycle, (ii) in the second one, named continuous simulations, the apparent plasticity at zero-stress is kept in memory for the next cycle. In other words, the continuous simulations consider the thermo-mechanical history seen by the material during the successive compressions. The internal intermolecular and network stresses are accumulated in different extents according to their dependence on the temperature and on the strain level reached. The results of the two kinds of simulations are presented in Figures 6 and7, respectively, and, with no surprise, the model is found to describe in a better manner the response by considering continuous simulations. The agreement between simulations and experiments is very good when the compression temperature increases. When the compression temperature decreases, the accumulation of the network stress stored in the material provokes an overestimation of the response for the lowest temperature. Even if the strain level was small enough to avoid orientation hardening, the relative contribution of the network stress increases with temperature. Xia and Tian (2019) showed that the strain hardening behavior of pre-deformed amorphous polymers is more accurately captured using a dissipative dynamics of molecular network orientation than using a temperature-dependent hardening modulus. The difference observed in Figure 7b at room temperature is thus potentially due to the way (from the rheological point of view) that the temperature-dependent hardening modulus (characterizing the molecular network orientation process) interacts with the temperature-dependent viscosity of the strain-hardening (characterizing the molecular network relaxation process). The thermally-activated strain recovery is related to an abrupt change in the ability of the polymer chain segments to rearrange towards equilibrium when the temperature traverses the glass transition. At a heating temperature of 433 K, Figure 8 provides the zero-stress creep response of polycarbonate specimens previously deformed at a strain of 70% and a temperature of 298 K. The recovery response is characterized by an initial rapid transient followed by the trend to reach a stabilization for which there is no significant strain change. It can be observed that the recovery response depends on the previous mechanical history in terms of strain-rate. The initial rapid transient and the equilibrium exhibit also a dependence on the heating temperature as shown in Figure 9. The experiments in Figures 8 and9 were replicated by the model using the constants calibrated in Section 3.2. The simulations were performed by a compression at the glassy state (at 298 K) up to a strain of 70% followed by an unload until zero-stress. Then, while the zero-stress was maintained constant, the remaining of the simulation was a temperature jump to the glass transition followed by a progressive heating at a rate of 0.05 K s -1 until the holding temperature of heating. Figures 8 and9 depict the ability of the model to predict the recovery response. It appears that the level of agreement between the predicted and measured recovery response depends on the previous applied strain-rate. However, the general trends provided by the model may be considered as relatively satisfactory. Interestingly, the model is able to capture the dependence of the response on thermo-mechanical history and recovery conditions. The equilibrium reached at the end of the zero-stress creep response does not lead to a complete recovery of the apparent plasticity. This incomplete strain recovery is highly dependent on coupled effects of strainrate and heating temperature and is reproduced by the model thanks to the consideration of the intrinsic strain-rate/temperature-dependent response coupled to the recovery-assisted network stress above Tg. Indeed, in our approach, viscoplastic intermolecular deformation and viscous network deformation are both strain-rate dependent and thermally-activated during the different steps of the loading-unloading-recovery process (Figure 4). The level of applied deformation is also an important condition of the mechanical history by accumulation of the two deformations. Its effect on the degree to which the strain is recovered is seen in the predictions provided in Figure 10. The model shows an important strain-rate dependency of the strain-level effect on the equilibrium state.

Concluding remarks

As concluding remarks, we have presented a physically-based model able to reproduce the elastoviscoplastic-viscohyperelastic response of amorphous thermoplastics in connection to the underlying thermo-mechanical mechanisms for below and above glass transition. The model successfully reproduced a series of experimental stress-strain data on polycarbonate and the abrupt change in response to a temperature variation over the entire range of strainrates tested. The material kinetics with temperature was then used as inputs in the internal stress of network orientation/relaxation to drive the thermally-activated strain recovery process and the predictive capabilities of the model were verified under various conditions. The generality endowed by the mathematical structure of the constitutive theory makes it applicable to the thermo-mechanical response of all amorphous thermoplastics and provides a useful tool to predict the thermally-induced strain recovery. In forthcoming works, MD simulations will be employed to better understand the deformation/recovery process in polymer systems in the final aim to propose enhanced continuum-based models.

Appendix A. Kinematics

The polymer medium is regarded as a homogeneous continuum body. As a basic quantity, the deformation gradient = ∂ ∂ F x X maps a material point from its initial position X to the current position x once the body is loaded. The determinant of the deformation gradient F is noted ( )

det 0 J = > F and its time derivative is = & F LF in which = ∂ ∂ L v x is the spatial velocity gradient with t = ∂ ∂ v x .
The initial configuration of the continuum body is at a homogeneous absolute temperature 0 θ taken as reference and the current configuration is at a homogeneous absolute temperature θ different than 0 θ . Using the conceptual sequence of configurations as proposed by [START_REF] Lee | Elastic-plastic deformation at finite strains[END_REF] with the introduction of an intermediate stress-free thermal configuration at θ , the thermal deformation effects can be combined to the mechanical deformation of the continuum body and the deformation gradient F can take a multiplicative form. It is decomposed into a mechanical part F M related to the isothermal stress-producing contribution and a thermal part θ F related to the stress-free thermal expansion of the continuum body:

M θ = F F F (A1) in which M θ = ∂ ∂ F x X and θ θ = ∂ ∂ F X X
= = F R U V R I I I I I (A7) h h h h h = = F R U V R N N N N N and v v v v v = = F R U V R N N N N N (A8)
The intermolecular and network spatial velocity gradients, I

(A9) -1 -1 -1 h v h h h v v h = + L L L F F F F F F & & 123 14243 N N N N N N N N N (A10)
The inelastic spatial velocity gradients, vp L I and v L N , may be both decomposed as the sum of symmetric and skew-symmetric parts:

= = + L F F F F D W & (A11) -1 -1 v h v v h v v N N N N N N N = = + L F F F F D W & (A12
) in which D is the stretching rate tensor (symmetric part) and W is the spin tensor (skewsymmetric part):

( )

1 2 T = + D L L and
( )

1 2 T = - W L L (A13)
It is assumed that the inelastic flow is volume preserving, i.e. Figure 2. Molecular events in an amorphous network loaded in compression followed by the thermally-induced shape-memory cycle (the internal stress includes intermolecular and network stresses) and one-dimensional visualization of the model; the amorphous polymer deformation F is split into a thermal part θ F related to the stress-free thermal expansion and a mechanical part F M related to the isothermal stress-producing contribution consisting in an elastoviscoplastic intermolecular branch in parallel with a viscohyperelastic network branch.

  average number of chains per unit volume) and the average number of rigid links rl N between physical molecular chain entanglements. As shown in Figure3c, the variation with temperature of the two parameters differs on both sides of g T . According to experimental observations of G'sell and Souami (1997) on the hardening modulus variation with temperature of amorphous networks, the parameter r C decreases below g T and increases above g T , whereas the parameter rl N increases below g T and is constant above g T . The following kinetics is then retained[START_REF] Richeton | Modeling and validation of the large deformation inelastic response of amorphous polymers over a wide range of temperatures and strain rates[END_REF]: , a and b are constants in the kinetics (17) and (18). The different behavior of r C and rl N on both sides of the glass transition is due to relative importance of internal energy and entropic effect in the glassy state and in the rubbery state.[START_REF] Richeton | Modeling and validation of the large deformation inelastic response of amorphous polymers over a wide range of temperatures and strain rates[END_REF] used the same kinetics (

F

  The intermolecular Cauchy stress σ I is expressed as a function of the corresponding strain by the following constitutive relationship with linear spring analog: represents the elastic volume change and ( ) e ln V I is the Hencky elastic strain (see Appendix A). The elastic stiffness tensor e C I is expressed in terms of the Young's modulus E and the Poisson's ratio ν as follows:

NFFN

  The deformation h in the hyperelastic spring drives Eq. (24) via the interaction with the deformation v N in the viscous damper (see Appendix A).The viscoplastic stretching rate tensor vp I D , and its alter ego, the viscous stretching rate tensor v N D , are given by the following flow rules: γ& are the accumulated viscoplastic and viscous strain rates, I σ and N σ are the effective shear stresses related to the deviatoric part of the intermolecular Cauchy stress and of the network Cauchy stress, respectively:

  By considering volume preserving and irrotationality of the inelastic flow (see Appendix A), the intermolecular stress I σ may be obtained by combining Eqs. (

  are defined into the two successive arrangements, θ X being the position of the material point in the stress-free thermal configuration. The determinant of the deformation gradient is M J J J θ = in which M J and J θ represent the intrinsic mechanical and thermal volume changes of the continuum body: 14243 (A3) in which ML is the isothermal stress-producing part and θ L is the stress-free thermal part.Due to the parallelism of the two resistances (see Figure2), the intermolecular and network deformation gradients, F I and F N , are both equal to the mechanical deformation gradient F M : relaxed configuration at θ , the two deformation gradients, F I and F N , can be in turn multiplicatively decomposed into elastic (or hyperelastic) and inelastic (viscoplastic or viscous) parts according to Lee decomposition: can be further decomposed into rotation R and stretching U or V (right or left, respectively) components:

FF

  frequently used hypothesis in the kinematics is irrotationality of the inelastic flow[START_REF] Gurtin | The decomposition F= FeFp, material symmetry, and plastic irrotationality for solids that are isotropic-viscoplastic or amorphous[END_REF], i.e. , of the intermolecular and network deformation gradients are then extracted from: , are then extracted from Eqs. (A5) and (A6). The tensors vp D I and v D N defining the specificity of the model, with its strain-rate and temperature dependence, are constitutively outlined in the Subsection 2.2.

Figure 1 .

 1 Figure 1. Thermally-induced shape-memory cycle of an amorphous polymer: deformation below g T ( I II → : loading, II III → : unloading), natural cooling at 0 T ( III IV → ) and recovery of shape and remanent strain upon heating above g T ( IV V → ).

  Intermolecular and network properties of polycarbonate as a function of temperature: (a) Young's modulus (1: 0.01 s -1 , 2: 0.001 s -1 , 3: 0.0001 s -1 ), (b) yield stress (1: 0.01 s -1 , 2: 0.001 s -1 , 3: 0.0001 s -1 ), (c) hardening modulus, average number of rigid links and viscosity constant, (d) work hardening rate-stress response at 0.001 s -1 (1: 298 K, 2: 333 K, 3: 363 K, 4: 393 K). The dashed lines represent the model and the symbols (the solid lines in (d)) represent the experiments.

Figure 4 .Figure 7 .

 47 Figure 4. Loading-unloading-recovery process and corresponding viscoplastic intermolecular stretch as a function of viscous network stretch under different strain-rates, 1: 0.01 s -1 , 2: 0.001 s -1 , 3: 0.0001 s -1 , the steps I V →

The shape-memory behavior is widely studied by a sequence of heating, deformation (in the rubbery state), cooling (with the deformation maintained) and unloading (in the glassy state); the heating closes the shapememory loop. The latter shape-memory cycle procedure is called hot programming within the shape-memory community(Dai et al.

, 2020). By contrast, a cold programming[START_REF] Dai | Modeling the thermo-mechanical behavior and constrained recovery performance of cold-programmed amorphous shape-memory polymers[END_REF] is preferred here. It consists to a sequence of deformation, unloading (in the glassy state), cooling and heating (in the rubbery state under zerostress creep) in order to highlight the importance of the elastoviscoplastic behavior in the glassy state.

The following notation is used throughout the text. Tensors and vectors are denoted by normal boldfaced letters and italicized boldfaced letters, respectively, while scalars and individual components of vectors and tensors are denoted by normal italicized letters. The superposed dot expresses the time derivative. The term tr denotes trace. The superscript T denotes the transpose.

In the specific case of polycarbonate, the α and β processes are associated to the rotation of the polymer main chain and the phenyl group, respectively.

& & & & & (2)

& & & &(4)in which the terms

& & & & g

Typical values of polycarbonate are used.

As one of the most successful theory in polymer physics,[START_REF] Doi | The Theory of Polymer Dynamics[END_REF] introduced the concept of the reptational Brownian motion: A macromolecule can be thought as trapped in a tube-like region where it can move in a snake-like fashion while having its lateral movements restricted. In an amorphous thermoplastic, the latter is due to the presence of neighboring chains and physical entanglements in the macromolecular network.

The superscripts e, vp, h and v denote the elastic, viscoplastic, hyperelastic and viscous contributions, respectively.

The reference glass transition is 423 K and the reference melting temperature is 493 K. These are typical values for polycarbonate[START_REF] Neki | Morphology-property studies of amorphous polycarbonate[END_REF].

The thermal expansion coefficient, which may be deduced from the slope of the dilatation-temperature curve (Ovalle-Rodas et al., 2015), is ( ) µ θ = 70.2 × 10 -6 m/mK below the glass transition and ( ) µ θ = 70.2 × 10 -6m/mK + 0.113/3Tg above the glass transition.

& (A14)1


Reza Saied

Ahmad Naddaf
email: naddaf@pasteur.ac.ir.highlights

Ahmad Mahmoudi

Mahdi Ghasemi

Ali Rohani

Seyyed Payman Mohammadi

Ziapour

Hesam Amir

Ehsan Nemati

Mostafavi

Saied Reza Naddaf
email: saiedrezanaddaf@gmail.com

Ahmad Ghasemi

Mahdi Rohani

Ali Mohammadi

Payman Ziapour

Nemati

Ehsan Mostafavi

Naddaf

Infection of hard ticks in the Caspian Sea littoral of Iran with Lyme borreliosis and relapsing fever borreliae

Keywords: Lyme borreliosis, Relapsing fever, Hard ticks

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

Introduction

The genus Borrelia in the family Spirochaetaceae comprises two distinct groups, each containing multiple arthropod-borne species that infect vertebrates. One group includes the species that cause relapsing fever (RF) and the other encompasses the causative agents of Lyme borreliosis (LB), collectively known as Borrelia burgdorferi sensu lato (s.l.) complex [START_REF] Parola | Ticks and tickborne bacterial diseases in humans: an emerging infectious threat[END_REF]. Except for the louse-adapted Borrelia recurrentis, the majority of relapsing fever borreliae (RFB) pathogenic to humans are vectored by the argasid soft ticks [START_REF] Talagrand-Reboul | Relapsing fevers: neglected tick-borne diseases[END_REF][START_REF] Cutler | Emerging borreliae -Expanding beyond Lyme borreliosis[END_REF]. Additionally, Ixodes and other hard ticks maintain Borrelia species such as Borrelia miyamotoi, Borrelia lonestari, Borrelia theileri, and some yet undescribed species that genetically group with RFB [START_REF] Talagrand-Reboul | Relapsing fevers: neglected tick-borne diseases[END_REF][START_REF] Furuno | Epidemiological study of relapsing fever borreliae detected in Haemaphysalis ticks and wild animals in the western part of Japan[END_REF][START_REF] Cutler | Emerging borreliae -Expanding beyond Lyme borreliosis[END_REF].

Soft tick-borne relapsing fever (STBRF) is endemic to many areas in Iran [START_REF] Karimi | Borrelia persica and B. baltazardi sp. nov.: experimental pathogenicity for some animals and comparison of the ultrastructure[END_REF][START_REF] Karimi | Relapsing fever and its epidemiology (in Persian)[END_REF][START_REF] Masoumi Asl | The epidemiology of tick-borne relapsing fever in Iran during 1997-2006[END_REF][START_REF] Naddaf | Tickborne relapsing fever in southern Iran, 2011-2013[END_REF]. Animal pathogenicity and molecular characterization assays indicated various RF borreliae species in Iranian Ornithodoros soft ticks [START_REF] Karimi | Relapsing fever and its epidemiology (in Persian)[END_REF][START_REF] Masoumi Asl | The epidemiology of tick-borne relapsing fever in Iran during 1997-2006[END_REF][START_REF] Naddaf | Phylogenetic analysis of the spirochete Borrelia microti, a potential agent of relapsing fever in Iran[END_REF][START_REF] Rafinejad | Detection of Borrelia persica infection in Ornithodoros tholozani using PCR targeting rrs gene and Xenodiagnosis[END_REF].

Borrelia persica, transmitted by Ornithodoros tholozani, is the primary cause of STBRF in the country [START_REF] Karimi | Borrelia persica and B. baltazardi sp. nov.: experimental pathogenicity for some animals and comparison of the ultrastructure[END_REF][START_REF] Karimi | Relapsing fever and its epidemiology (in Persian)[END_REF][START_REF] Shirani | A case of canine borreliosis in Iran caused by Borrelia persica[END_REF]. In regions where O. tholozani does not occur, epidemiological and molecular evidence suggests Borrelia microti [START_REF] Janbakhsh | The nature of sporadic cases of relapsing fever in Kazeroun area, southern Iran[END_REF][START_REF] Masoumi Asl | The epidemiology of tick-borne relapsing fever in Iran during 1997-2006[END_REF][START_REF] Naddaf | Phylogenetic analysis of the spirochete Borrelia microti, a potential agent of relapsing fever in Iran[END_REF], and closely related species, the presumed ecotypes of Borrelia duttonii, as the cause of the disease [START_REF] Naddaf | Tickborne relapsing fever in southern Iran, 2011-2013[END_REF][START_REF] Naddaf | Relapsing fever causative agent in Southern Iran is a closely related species to East African borreliae[END_REF].

At present, not much data on LB and the relevant spirochetes are available from Iran. LB is the most prevalent tick-borne infection in the temperate regions of North America, Europe, and some parts of Asia [START_REF] Stanek | Lyme borreliosis[END_REF]. The causative agents are transmitted among animals and to humans via the infective bites of hard ticks belonging to the genus Ixodes (Acari: Ixodidae) [START_REF] Brites-Neto | Tick-borne infections in human and animal population worldwide[END_REF]. Unlike RFB, no evidence of transovarial transmission of B. burgdorferi s.l. is available, and they seem to perpetuate in cycles involving vertebrate reservoirs and tick vectors [START_REF] Breuner | Lack of evidence for transovarial transmission of the Lyme disease spirochete Borrelia mayonii by infected female Ixodes scapularis (Acari: Ixodidae) Ticks[END_REF]. White-footed mice (Peromyscus leucopus) are the principal reservoir host for B. burgdorferi sensu stricto (s.s.), the primary etiological agent of LB in eastern North America [START_REF] Roy-Dufresne | Poleward expansion of the white-footed mouse (Peromyscus leucopus) under climate change: implications for the spread of lyme disease[END_REF]. Whereas, in Eurasia, a variety of small mammals, including field mice (Apodemus spp.), the bank vole (Myodes glareolus), shrews (Sorex spp.), red squirrels (Sciurus vulgaris), and several bird species serve as reservoir hosts for variants of B. burgdorferi s.l. [START_REF] Matuschka | Capacity of European animals as reservoir hosts for the Lyme disease spirochete[END_REF][START_REF] Gern | European reservoir hosts of Borrelia burgdorferi sensu lato[END_REF][START_REF] Mysterud | Genospecies of Borrelia burgdorferi sensu lato detected in 16 mammal species and questing ticks from northern Europe[END_REF].

The Hyrcanian mixed forest of Caspian Sea littoral, Iran, is the southernmost distribution range of Ixodes ricinus ticks in Eurasia [START_REF] Stanek | Lyme borreliosis[END_REF]. Many studies have reported this species parasitizing livestock in this area [START_REF] Sofizadeh | Hard tick species of livestock and their bioecology in Golestan province, North of Iran[END_REF][START_REF] Razmi | Prevalence of ixodid ticks on cattle in Mazandaran province, Iran[END_REF]. Moreover, the possible rodent reservoirs of the LB agents, e.g., the members of the genus Apodemus, like A. hyrcanicus, A. uralensis, and A. witherbyi prevail in the area [START_REF] Ghorbani | Rodent fauna of the western Golestan Province in northeast Iran[END_REF][START_REF] Esfandiari | Epidemiological distribution of rodents as potent reservoirs for infectious diseases in the provinces of Mazandaran, Gilan and Golestan, Northern Iran[END_REF].

Despite the occurrence of I. ricinus and many other hard tick species in the Caspian Sea littoral of Iran, no data on the borrelial infection of these ticks is available from this area. Here, for the first time, we report the infection of the I. ricinus and other hard ticks in the Caspian Sea littoral, Iran, with LB and RF borreliae. We also investigated the presence of borrelial DNA in the common rodents and small mammals, the potential reservoirs of tick-borne spirochetes.

Materials and Methods

Study area

The study area included different geological landscapes, e.g., plains, grasslands, foothills, mountains, river-banks, and forests in Golestan, Mazandaran, and Gilan provinces, Caspian Sea littoral, Iran. No region in Iran receives as much precipitation as this area. According to the Köppen-Geiger classification, of the five main climate groups, three groups, B (dry), C (temperate), and D (continental), occur in Iran. The Caspian Sea littoral is within the C region with subdivisions Csb (warm-summer Mediterranean), and Csa (hot-summer Mediterranean climate), in the Golestan and Mazandaran, and a combination of these two and Cfa (humid subtropical climate) in Gilan [START_REF] Kottek | World Map of the Köppen-Geiger climate classification updated[END_REF][START_REF] Raziei | Köppen-Geiger climate classification of Iran and investigating its changes during 20th century[END_REF]. These provinces encompass parts of Hyrcanian mixed forests, a green belt of predominantly temperate deciduous forests stretching over the northern slopes of the Alborz mountain range, down to the southern borders of the Caspian Sea. These forests are about 20-50 million years old and serve as refugia for many Arcto-Tertiary relict species [START_REF] Zohary | Geobotanical foundations of the Middle East[END_REF][START_REF] Leestmans | Le refuge caspiens et son importance en biogéographie[END_REF] (Fig. 1).

Tick collection

The ticks were removed by fine-pointed, stainless-steel tweezers from mammalian hosts, including sheep, goats, cattle, camels, horses, dogs, donkeys, rodents, and hedgehogs. For the collection of the unfed ticks (larvae, nymphs, and adults), we used the dragging method in the grasslands at the border of forested areas. The collected ticks in falcon tubes, along with the background information, including the host, locality name, date of collection, and the collector's name, were transferred to the Parasitology Laboratory in Amol Research Center, Pasteur Institute of Iran. The specimens were transferred to 70% alcohol, separated based on gender, and identified morphologically using taxonomic keys [START_REF] Walker | The genus Rhipicephalus (Acari, Ixodidae): a guide to the brown ticks of the world[END_REF][START_REF] Walker | Ticks of domestic animals in Africa. A guide to identification of species[END_REF][START_REF] Apanaskevich | The genus Hyalomma Koch, 1844. ii. Taxonomic status of H. (Euhyalomma) anatolicum Koch, 1844 and H. (E.) excavatum Koch, 1844 (Acari, Ixodidae) with redescriptions of all stages[END_REF][START_REF] Apanaskevich | The genus Hyalomma Koch, 1844: V. Re-evaluation of the taxonomic rank of taxa comprising the H. (Euhyalomma) marginatum Koch complex of species (Acari: Ixodidae) with redescription of all parasitic stages and notes on biology[END_REF][START_REF] Apanaskevich | The genus Hyalomma Koch, 1844. X. redescription of all parasitic stages of H. (Euhyalomma) scupense Schulze, 1919 (= H. detritum Schulze) (Acari: Ixodidae) and notes on its biology[END_REF][START_REF] Hosseini-Chegeni | The updated list of ticks (Acari: Ixodidae & Argasidae) occurring in Iran with a key to the identification of species[END_REF][START_REF] Hosseini-Chegeni | The Iranian Hyalomma, with a key to the identification of male species[END_REF]. The specimens were kept at 4°C until used for DNA extraction.

Rodent trapping

The live traps were set in different localities before dusk and collected the next day early morning. The selection of localities was based on the area's topography to cover the habitats of different species. In grasslands, plains, meadows, and farms, the traps were placed at the entrance of active burrows, while for the forest-dwelling species, the traps were randomly placed deep into the wooded areas or along the riverbanks. The captured animals were examined on-site for ectoparasites and then transferred to the laboratories affiliated to one of the state universities in the three provinces. The animals were euthanized according to the ethical codes of the Pasteur Institute of Iran (Reference No. IR.PII.REC.1395.29), and blood samples were collected from the heart of the animals. The animal sera and pieces of the spleen, liver, lung, and muscles were kept in alcohol 96% at -20° C for future investigations. The specimens were compared with the vouchers deposited at the Mohammad Hanifi's Health Museum (MHHM), Akanloo Research Center, Pasteur Institute of Iran, and identified based on the morphological keys available in the literature [START_REF] Kryštufek | Mammals of Turkey and Cyprus: Rodentia I: Sciuridae, Dipodidae, Gliridae, Arvicolinae[END_REF][START_REF] Ghorbani | Rodent fauna of the western Golestan Province in northeast Iran[END_REF][START_REF] Darvish | Faunistic and taxonomic study of Rodents from northwestern Iran[END_REF].

DNA extraction

Individual and pooled ticks or 25-50 mg of spleen and liver of the rodents were ground to powder in 2 ml tubes containing liquid nitrogen by using a pestle grinder. Ticks of the same species, sex, and developmental stage (nymph, and adult) were pooled in groups of 3-22 according to the host, location/date of collection, and the size. DNA extraction from the animals' visceral tissues was performed using a commercial kit (QIAGEN, Hilden, Germany) according to the manufacturer's instructions, and from ticks with a potassium acetate procedure successfully practiced by others [START_REF] Rodríguez | An alternative and rapid method for the extraction of nucleic acids from ixodid ticks by potassium acetate procedure[END_REF] with minor modifications. In the latter procedure, I) 500 µL of lysis buffer [0.1 M Tris-HCl (pH 8.25), 0.05 M EDTA, 0.2 M sucrose, 0.5% SDS] supplemented with proteinase K (10 mg/mL) was added to the ground specimens. II) the suspensions were incubated at 56 o C for 24 h. III) following the addition of 120 µl of 5M potassium acetate and incubation on ice for 10 min, the suspensions were centrifuged at 12 000 ×g for 10 min to save supernatants. IV) 35 µl of 4M sodium acetate and 1000 μl ethanol (96-100%) were added to the recovered supernatants followed by incubation on ice for 10 min and then centrifugation at 12 000 ×g for 20 min. V) the supernatants were decanted, and the pellets were gently washed with 500 μl 70% ethanol to remove residual salts. VI) the remaining alcohol was allowed to evaporate at room temperature for 24 h, and then 200 μl double distilled water (DDW) was added to the tubes and stored at -20°C until used.

qPCR

The whole genomic DNAs extracted from ticks and the animal's viscera were screened for borrelial DNA by targeting a 136-bp sequence using the primers and the probe (Table 1) specific to the 16S rRNA region of the Borrelia genus [START_REF] Ornstein | A reverse transcriptase-polymerase chain reaction assay of Borrelia burgdorferi 16S rRNA for highly sensitive quantification of pathogen load in a vector[END_REF]. The final 20 μl reactions contained 10 μl master mix 2X (Amplicon, Denmark), 500 nM of each primer, 200 nM probe, and 4 μl (50 nM) of the template DNA. Amplifications were performed in a Rotor-Gene 6000 instrument (Corbett Life Science, Sydney, Australia) programed for an initial denaturation at 95 °C for 10 min, followed by 45 cycles of denaturation at 95 °C for 15 s, and annealing/extension at 60 °C for 60 s. DDW as no template control (NTC) was included in all assays to ensure that the reagents were free of DNA contamination. Ten-fold serial dilutions of B. burgdorferi s.s. DNA (Amplirun® Borrelia DNA control), containing approximately 7200, 720, 72, 7.2, 0.72, and 0.072 copy numbers per reaction of bacterial DNA, was used in the control reactions to measure the limit of detection.

nested PCR and sequencing

The positive qPCR specimens were further analyzed by nested PCRs (nPCR) using two sets of primers complementary to flaB and the non-coding rrs-rrlA spacer (Table 1). The 25 µl reactions contained 12.5 µl 2X master mix (Ampliqon, Denmark), 10 pm of each of forward and reverse primer, 200 µM of dNTPs, 2 µl of template DNA, and DDW to the final volume. The first amplification program for both target fragments began with an initial denaturation at 95°C for 15 min followed by 40 cycles of 30 s at 95°C, 45 s at 50°C, and 60 s at 72°C and a final extension of 7 min at 72°C. The PCR products, diluted in a ratio of 1/20 with DDW, were used in the second amplification step with the same program except for increasing the extension temperatures to 54°C for flaB and 55°C for rrs-rrlA. The PCR products were resolved on 1% agarose gel, stained with DNA green viewer (Parstous Biotechnology, Mashhad, Iran) and visualized under UV. The amplicons were sequenced in both directions with the same primers used in the second amplification by the Sanger method in a commercial Company (Genomin, Tehran, Iran).

BLAST analysis

Ambiguous sites of the generated sequences (≈600 bp) were identified and resolved using CodonCode Aligner software (CodonCode Corp.). The sequences were compared with similar sequences of Borrelia species available in the GenBank database by using the Basic Local Alignment Search Tool (BLAST) analysis (http://blast.ncbi.nlm.nih.gov), and similarities were obtained.

Phylogenetic analysis

Overall, 41 sequences, including seven from the present study and 34 from the GenBank database (Table 2), were aligned with the Clustal W algorithm [START_REF] Thompson | CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice[END_REF], as implemented in BioEdit software ver.7.0.9 [START_REF] Hall | BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT[END_REF]. The best-fitting model of nucleotide substitution was separately determined for the Maximum Likelihood (ML) and Bayesian Inference (BI) analyses based on Akaike (AIC = TIM3+G) and Bayesian (BIC = TPM3uf+G) Information Criterion, respectively using jModeltest program [START_REF] Posada | jModelTest: phylogenetic model averaging[END_REF]. The ML analysis was performed to construct the molecular tree using RaxML v. 8.2.10 (Stamatakis, 2014) software with the bootstrap probability (BP) values obtained from 1000 replicates. The BI tree was constructed using MrBayes v3.2 [START_REF] Ronquist | MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space[END_REF]. Two independent analyses were performed with a Markov chain Monte Carlo (MCMC) model, which started from random trees and continued for 4×10 6 generations. The trees were sampled every 100 th generation after the first 25% of trees were discarded as burn-in, and Bayesian posterior probabilities (BPP) was used to assess branch support of the BI tree. BP ≥80 and BPP ≥ 0.95 were considered as "good" supports. The generated trees were rooted using the midpoint approach and visualized using FigTree v.1.4.2.

GenBank submission

The sequences generated in this study were deposited in the GenBank under the accession numbers MN958341 to MN958351.

Results

Tick species

We collected 1094 and 1475 ticks from Mazandaran in spring and autumn, respectively, and 1569 from Golestan, and 176 from Gilan. Overall, 13 different species belonging to four genera were identified in the study area. In Mazandaran, we found ten species with I. ricinus as the most collected species (71.2%), followed by Rhipicephalus sanguineus sensu lato (s.l.) (12.3%), and Rhipicephalus annulatus (7.5%). In Golestan, we identified nine species, Rhipicephalus turanicus was the primary species (64.7%) followed by Rh. sanguineus s.l. (21.5%), and I. ricinus (5.1%). Most of the ticks were removed from livestock and a few from hedgehogs. In Gilan Province, due to the application of acaricides, we collected much fewer specimens, mostly I. ricinus (76.1%). In Mazandaran, we succeeded in collecting some nymphs of I. ricinus and Rh.

annulatus ticks from livestock. However, despite several attempts, we could not collect unfed ticks by the dragging method. Table 3 shows the details of the collected ticks in the study area.

Rodent species

In three provinces, Golestan, Mazandaran, and Gilan, we captured a total of 39 small mammals comprising six and three species from the orders Rodentia (n=33), and Eulipotyphla (n=6) (Table 4 ).

qPCR

Out of 509 DNA specimens of the individual and pooled ticks, 71 belonging to I. ricinus ticks or the members of the Rhipicephalus genus (Rh. sanguineus s.l., Rh. turanicus, Rh. bursa, Rh. annulatus) became positive for Borrelia DNA by qPCR with the threshold (Ct ) value ≤ 40 (Table 5). The mean Ct values ± SD for the reactions containing ≈7200, ≈720, ≈72, ≈7.2, and ≈0.72, copy numbers were 29.3±0.8, 32.6±0.8, 36.6±0.9, 40.5±0.4, and 42.2±0.1, respectively; no measurable Ct value was obtained for the NTCs and the reactions containing ≈0.072 copy number. Other tick species and the animals were negative for the borrelial DNA.

nested PCR

Among 71 qPCR positive specimens, amplification of flaB yielded a ≈600-bp band in 19 specimens belonging to I. ricinus (n=8), Rh. sanguineus s.l. (n=8), Rh. turanicus (n=2) and Rh. bursa (n=1). Also, in 67 qPCR positive specimens, PCR amplification of rrs-rrlA spacer revealed a ≈530 bp in 42 specimens from Rh. sanguineus s.l. (n=17), Rh. turanicus (n=15), Rh. bursa (n=1), and I. ricinus (n=9) ticks, and a ≈1000 bp amplicon in a pooled I. ricinus DNA (Table 5).

BLAST analysis

The flaB sequences generated from the individual or pooled I. ricinus ticks showed the highest identity (99.50%-100%) with the sequences representing Borrelia bavariensis PBi (acc. No. 6). The ≈530-bp rrs-rrlA sequences from Rhipicephalus spp. and I. ricinus ticks matched only ≤ 40% of subject sequences in the GenBank database and did not resolve the identity of Borrelia species. The only rrs-rrlA 940-bp sequence from a pooled I. ricinus in Golestan matched 98.19% with B. valaisiana Tom4006 (Acc. No.CP009117) [START_REF] Kurilshikov | Complete Genome Sequencing of Borrelia valaisiana and Borrelia afzelii Isolated from Ixodes persulcatus Ticks in Western Siberia[END_REF].

Phylogenetic analysis

The BI and ML trees were completely congruent; hence, only the ML tree is shown here. The

2).

Discussion

The Caspian seal littoral has a rich fauna of hard ticks; of the 36 ixodid tick species identified in Iran [START_REF] Hosseini-Chegeni | The updated list of ticks (Acari: Ixodidae & Argasidae) occurring in Iran with a key to the identification of species[END_REF], at least 19 species, including I. ricinus, the notorious vector of LB in Eurasia [START_REF] Razmi | Prevalence of ixodid ticks on cattle in Mazandaran province, Iran[END_REF][START_REF] Sofizadeh | Hard tick species of livestock and their bioecology in Golestan province, North of Iran[END_REF], occur in this area.

In this study, by analyzing the flaB sequence, we, for the first time, report the infection of I. 1) in a 2% agarose gel [START_REF] Wodecka | A comparative analysis of molecular markers for the detection and identification of Borrelia spirochaetes in Ixodes ricinus[END_REF]. In our study, the same primers yielded only ≈530 bp and ≈1000 bp amplicons. The ≈530-bp sequences, despite transparency, did not resolve much of the identity of Borrelia species and the only 940-bp sequence generated from pooled I. ricinus ticks (Table 5.

No. 55) matched 98.19% with B. valaisiana Tom4006 (Acc. No.CP009117) (Fig. 1).

qPCR detected borrelial DNA in 71 out of 501 (≈14%) tick specimens; however, only 19 and 43 specimens became positive by PCR assays that amplified the flaB, and rrs-rrlA spacer, respectively. Contrary to our findings, the primers complementary to the flaB showed a much higher sensitivity than primers complementary to the rrs-rrlA spacer in detecting borrelial DNA in Ixodes ticks [START_REF] Wodecka | A comparative analysis of molecular markers for the detection and identification of Borrelia spirochaetes in Ixodes ricinus[END_REF]. The difference in the positivity rate of the qPCR and nPCR assays in our study might be due to the higher sensitivity of the former one, which previously detected DNA concentrations equivalent to as low as one-thousandth of a single B. burgdorferi cell [START_REF] Ornstein | A reverse transcriptase-polymerase chain reaction assay of Borrelia burgdorferi 16S rRNA for highly sensitive quantification of pathogen load in a vector[END_REF]. Similar studies have also reported such discrepancies in the sensitivity of qPCR and nPCR, mainly when different loci are targeted [START_REF] Gil | Identification of a new Borrelia species among small mammals in areas of northern Spain where Lyme disease is endemic[END_REF]. Here, our qPCR could detect as low as ≈0.72 copy number in reactions. None of the 37 specimens that became negative by qPCR produced an amplicon in the PCR assays.

The generated flaB sequences exhibited 99.50%-100% identities with European LB agents representing B. afzelii, B. garinii, B. valaisiana, or B. bavariensis (Table 6). In Europe, B. afzelii and B. garinii accounts for the majority of LB cases, whereas, in Asia, B. garinii is the primary cause of the disease [START_REF] Stanek | Lyme borreliosis[END_REF]. Borrelia afzelii mostly remains localized in the skin tissue, while B. garinii, in most cases, is neurotropic [START_REF] Margos | Population genetics, taxonomy, phylogeny and evolution of Borrelia burgdorferi sensu lato[END_REF][START_REF] Stanek | Lyme borreliosis[END_REF].

Other LB spirochetes, Borrelia bissettiae, Borrelia lusitaniae, Borrelia spielmanii, and B. valaisiana, are also known to cause human infections in Europe but are not significant pathogens [START_REF] Diza | Borrelia valaisiana in cerebrospinal fluid[END_REF][START_REF] Maraspin | Lyme borreliosis and Borrelia spielmanii[END_REF][START_REF] Collares-Pereira | First isolation of Borrelia lusitaniae 11 from a human patient[END_REF][START_REF] Rudenko | Molecular detection of Borrelia bissettii DNA in serum samples from patients in the Czech Republic with suspected borreliosis[END_REF].

Borrelia bavariensis, previously known as B. garinii OspA serotype 4, has a wide distribution in Eurasia [START_REF] Margos | Borrelia bavariensis sp. nov. is widely distributed in Europe and Asia[END_REF]. In Europe, B. bavariensis strains are the primary cause of Lyme neuroborreliosis [START_REF] Fingerle | Epidemiological aspects and molecular characterization of Borrelia burgdorferi s.l. from southern Germany with special respect to the new species Borrelia spielmanii sp[END_REF], and in Japan, most human infections are linked to this species [START_REF] Takano | Multilocus sequence typing implicates rodents as the main reservoir host of human-pathogenic Borrelia garinii in Japan[END_REF]. Ixodes persulcatus and I. ricinus transmit B. bavariensis in Asia and Europe, respectively [START_REF] Margos | Borrelia bavariensis: Vector switch, niche invasion, and geographical spread of a tick-borne bacterial parasite[END_REF]. It was hypothesized that range expansion of this species from Asia to Europe coincided with a change in the spirochete's tick vector, switching from I. persulcatus in Asia to I. ricinus in Europe [START_REF] Margos | Borrelia bavariensis: Vector switch, niche invasion, and geographical spread of a tick-borne bacterial parasite[END_REF]. Hence, the presence of this species in I. ricinus ticks of northern Iran suggests that the B. bavariensis strain found in this study is of European origin. Our flaB phylogeny provided further evidence for the range expansion of European hard-tick borreliae into Iran, e.g., the B. myiamotoi strains from north of Iran grouped closer to the European B. miyamotoi (100/1.00) than the American and Chinese strains (Fig. 2).

No official report of human LB infection in Iran is available; however, a first-tier ELISA with subsequent immunoblotting detected anti-B. burgdorferi antibodies in dogs in the Caspian Sea littoral [START_REF] Hanifeh | Seroprevalence, geographic distribution and risk factor analysis of Borrelia burgdorferi sensu lato in naturally exposed dogs of Iran[END_REF][START_REF] Hanifeh | Effect of climatic factors on canine lime borreliosis in Iranian Caspian sea littoral provinces[END_REF]. Turkey, in the northwest, is the nearest country to Iran, in which infection of I. ricinus ticks with B. burgdorferi s.l. was reported (Guner et al., 2003;[START_REF] Gunes | The prevalence of Borrelia burgdorferi sensu lato in Ixodes ricinus ticks in the Sinop region of Turkey[END_REF], and the LB in humans were identified [START_REF] Onal | Systematic review of Lyme disease in Turkey[END_REF].

In our study, borrelial DNA was detected only in the adult ticks, mostly parasitizing the sheep and cattle. The lack of Borrelia infection in I. ricinus nymphs might be due to inadequate numbers of the specimens. Also, the presence of various species of lizards in the area [START_REF] Smid | Annotated checklist and distribution of the lizards of Iran[END_REF], as reservoir-incompetent hosts, might play a cleansing role and reduce the density of infected I. ricinus nymphs [START_REF] Swei | Impact of the experimental removal of lizards on Lyme disease risk[END_REF]. Despite several attempts, we were unable to catch the immature I. ricinus specimen by dragging in 2016 and 2017. However, we succeeded in collecting a limited number of nymphs in the 2018 autumn (unpublished data).

We did not detect LB agents in the viscera of the rodents and other small mammals. Borrelia DNA, following dissemination of the LB infection in the vertebrate host, generally become detectable in various tissues, including the viscera [START_REF] Sertour | Infection kinetics and tropism of Borrelia burgdorferi sensu lato in mouse after natural (via ticks) or artificial (needle) infection depends on the bacterial strain[END_REF], and viable spirochetes can be cultured from these tissues. In naturally and experimentally-infected P. leucopus mice and other small mammals, visceral tissues, including the spleen and liver, proved reliable for detecting LB spirochetes [START_REF] Anderson | Seasonal prevalence of Borrelia burgdorferi in natural populations of white-footed mice, Peromyscus leucopus[END_REF][START_REF] Barthold | Experimental Lyme arthritis in rats infected with Borrelia burgdorferi[END_REF]; however, they were not as much rewarding as the ear tissues [START_REF] Sinsky | Ear punch biopsy method for detection and isolation of Borrelia burgdorferi from rodents[END_REF]. Hence, in our study, by analyzing only the animals' visceral tissues, we might have missed some positive animals.

The RF-like sequences from Rhipicephalus ticks were identical and exhibited the highest identity (98.97%) with the Borrelia sp. HM (acc. No. LC170031) infecting Ha. megaspinosa ticks and wild boar in Japan [START_REF] Furuno | Epidemiological study of relapsing fever borreliae detected in Haemaphysalis ticks and wild animals in the western part of Japan[END_REF] (Table 6). However, in the phylogeny, the representative sequence (Borrelia sp. Ir-Maz173) emerged as a sister taxon to the Borrelia sp. identified in Rh. sanguineus s.l. of Portugal. (acc. No. KR677090) (Fig. 2).

Our phylogeny revealed two deeply divergent groups within the genus Borrelia, one, defined as the LB clade, encompassed all members of B. burgdorferi s.l. (100/1.00), and the other included the agents that clustered either in REP or RF clades (100/1.00). The controversial REP borreliae, B. turcica, and B. tachyglossi, grouped as a highly supported distinct clade (98/1.00), yet, shared a distant common ancestor with the RF borreliae (Fig. 2).

Further investigations with additional tick specimens, including the immature stages and employment of more resolving tools like multilocus sequence typing [START_REF] Hanincova | Multilocus sequence typing of Borrelia burgdorferi suggests existence of lineages with differential pathogenic properties in humans[END_REF][START_REF] Wang | Molecular typing of Borrelia burgdorferi[END_REF], are required to reveal other possible Borrelia species and population structure of the prevalent species in the area. Also, xenodiagnosis, i.e., deploying the natural arthropod vector for detection and recovery of the infectious agent [START_REF] Huegli | Apodemus species mice are reservoir hosts of Borrelia garinii OspA serotype 4 in Switzerland[END_REF], and using ear punch biopsy method [START_REF] Sinsky | Ear punch biopsy method for detection and isolation of Borrelia burgdorferi from rodents[END_REF] might identify the LB reservoirs with higher sensitivity.

The Caspian Sea littoral is endemic for other febrile infectious diseases like leptospirosis [START_REF] Zakeri | Molecular epidemiology of leptospirosis in northern Iran by nested polymerase chain reaction/restriction fragment length polymorphism and sequencing methods[END_REF][START_REF] Garshasbi | Leptospirosis in Caspian Sea littoral, Gilan Province, Iran[END_REF] that show a seasonal prevalence overlapping with tick activity and the associated tick-borne infections. Increased awareness of the possible LB and RF infections in this area can result in consideration of these agents in the differential diagnosis of infectious diseases, and the fever of unknown origin. [START_REF] Lescot | The genome of Borrelia recurrentis, the agent of deadly louse-borne relapsing fever, is a degraded subset of tick-borne Borrelia duttonii[END_REF]

Figure captions:

CP000013 ) ,

 , Borrelia garinii (acc. No. KR782234), Borrelia valaisiana (acc. No. MF150080), and Borrelia afzelii (Acc. No. MG944963), or the relapsing fever agent, B. miyamotoi (acc. No. KX6461990). All the sequences from the Rhipicephalus ticks were identical and exhibited the highest identity (98.97%) with an RF-like Borrelia (acc. No. LC170031) infecting Haemaphysalis megaspinosa ticks in Japan. They also shared 96.73% identity with an isolate in Rh. sanguineus of Portugal (acc. No. KT364297), and B. theileri ( acc. No. KF569936) of Rhipicephalus geigyi in Mali, and 96.04% with B. lonestari of Amblyomma americanum in USA (Table

  constructed phylogenetic trees based on the flaB sequences revealed two main groups, with high supports (BP=100, BPP=1.00), representing RF-like and LB borreliae. The RF-like borreliae showed a further relatively deep substructuring into two main clades, RF borreliae (RFB), and the recently defined reptile-associated (REP) borreliae. Here, REP borreliae, represented by Borrelia turcica (Acc. No. AB109245) and Borrelia tachyglossi (Acc. No. KY586966) appeared as the sister lineage of the RF borreliae with high levels of BP and BPP supports (100/1.00). The generated LB sequences grouped with either of the four species of the B. burgdorferi s.l. complex, including B. bavariensis (Acc. No. DQ650333), B. garinii IM710 (Acc. No. KR782234), B. afzelii clone 2-7 (Acc. No. KX646195), and B. valaisiana Bv-V4 (Acc. No. MF150080). The hard tick RF borreliae, including the sequences generated in the present study clustered separately from soft tick-associated RF species. The RF-like Borrelia detected in Rhipicephalus ticks grouped with two hard ticks borreliae from Portugal (Acc. No. KR677090) and Japan (Acc. No. LC170031) (85/1.00), closely related to B. theileri (KF569936), B. lonestari (AY850063), and another isolate from Portugal (KT364297). The RF-like borreliae identified in I. ricinus ticks grouped with a B. miyamotoi strain from Europe (Acc. No. KT948321), adjacent to the USA and China strains (Acc. Nos. CP017126 and KU749379) (Fig.

  ricinus ticks in the Caspian Sea littoral with LB spirochetes, B. bavariensis, B. garinii, B. afzelii, and B. valaisiana as well as the RF Borrelia, B. miyamotoi. The species B. bavariensis, B. valaisiana, and B. afzelii, along with B. miyamotoi, were detected in I. ricinus ticks of Mazandaran, mostly collected during the autumn, and B. garinii was identified in I. ricinus ticks of Gilan province. In Golestan, no flaB sequence from I. ricinus ticks was obtained. Previously, a PCR with primers complementary to the rrs-rrlA spacer identified B. afzelii, B. garinii, and B. burgdorferi s.s. based on the appearance of target band sizes (Table

Figure 1 .

 1 Figure 1. Map of the study area and the localities in the Caspian Sea littoral from where hard ticks and rodents were collected. The light green shadow shows the Hyrcanian mixed forests. The abbreviations stand for the identified species; Bg, B. garinii; Bb, B. bavariensis; Bv, B. valaisiana; Ba, B. afzelii; Bm, Borrelia miyamotoi; Bsp, Borrelia sp.

Figure 2 .

 2 Figure 2. The ML tree constructed based on flaB sequences of Borrelia species. The sequences generated in this study are in red. Numbers on the nodes refer to bootstrap probability (BP) and Bayesian posterior probabilities (BPP) values, respectively. The support values <70 are not shown here. The scale bar indicates 0.02-estimated nucleotide substitutions per site. † DQ650333 in the GenBank refer to B. garinii, but in the article (PUBMED No. 21841027) represents B. bavariensis. The accession numbers for B. garinii (Ir-Gil 2), B. bavariensis (Ir-Maz 6), B. valaisiana (Ir-Maz22), B. afzelii (Ir-Maz122), B. miyamotoi (Ir-Maz57), B. miyamotoi (Ir-Maz123) and Borrelia sp (Ir-M173) in the GenBank database are MN958341, MN958342, MN958343, MN958344, MN958345, MN958347, and MN958349, respectively.

Table 1 :

 1 The primers and the probe used for the amplification of 16S rRNA, rrs-rrlA, and flaB sequences.

	Bacteria Target	Primers/probe	Sequence (5' to 3')	Fragment	Reference
					size bp
	Borrelia	16S rRNA 16SBOR-Fw	GGTCAAGACTGACGCTGAGTCA	136	(Ornstein and
	spp.		16SBOR-Rev	GGCGGCCACTTAACACGTTAG		Barbour, 2006)
			16SBOR-P	Fam-
				TCTACGCTGTAAACGATGCACACTTGG
				TG-BHQ-1
	Borrelia	rrs-rrlA	Outer	P95 : 5′-GGTTAGAGCGCAGGTCTG-3′	1014	(Ranka et al.,
	spp.			Pa: 5′-GGTATGTTTAGTGAGGG-3′		2004, Wodecka
			Inner	P97 : 5′-GATGTTCAACTCATCCTGGTCCC-	B. afzelii 435	et al., 2010)
				3′	B. garinii 532
				Pb: 5′-CGTACTGGAAAGTGCGGCTG-3′	B. burgdorferi
					s.s 950
	Borrelia	flaB	Outer	132f: 5′-TGGTATGGGAGTTTCTGG-3′	774	(Wodecka,
	spp.			905r: 5′-TCTGTCATTGTAGCATCTTT-3′		2007)
			Inner	220f: 5′-CAGACAACAGAGGGAAAT-3′	604
				823r: 5′-TCAAGTCTATTTTGGAAAG
				CACC-3′

Table 2 :

 2 Details of flaB sequences used in phylogenetic analysis.

	Borrelia Species	Host/isolation	Location	GenBank	Reference
		source		Acc. No.
	B. persica str. T	O. tholozani	Iran	JF708953	Unpublished
	B. microti str. IR-1	O. erraticus	Iran	JF825472	(Naddaf et al., 2012)
	B. duttonii str. Ly	Human	Africa	CP000976	(Lescot et al., 2008)
	B. crocidurae str.	O. sonrai	Mauritania CP003426	(Elbir et al., 2012)
	Achema
	B. hispanica str.Sp3	O. erraticus	Spain	GU357612	(Toledo et al., 2010)
	B. recurrentis str. A1	Human	Ethiopia	CP000993

Acknowledgment

We are grateful to the Center for Communicable Disease Control (Ministry of Health), Golestan University of Medical Sciences, Mazandaran University of Medical Sciences, and the Gilan University of Medical Sciences for providing us the facilities during the study. We thank our

Funding resources

This study was funded by the Pasteur Institute of Iran (Grant No. 917).

enthusiastic colleagues Hamed Hanifi and Seyed Adel Hosseini, for their assistance in the sample collection. We also thank Mrs. Masoumeh Khoshyar for preparing the map of the sampled locations.

Conflict of interest

Saied Reza Naddaf and Ahmad Mahmoudi equally contributed to this work as the first author.

The authors declare there is no issue to be conceived as a conflict of interest related to this article. Ornithodoros; C., Calomys; I, Ixodes; A., Amblyomma; Rh., Rhipicephalus; Ha., Haemaphysalis; B., Bothriocroton; Hy., Hyalomma