Adria Roig

Sai Kishor Kothakota

Narcis Miguel

Pierre Fernbach

Enrico Mingo Hoffman

Luca Marchionni

On the Hardware Design and Control Architecture of the Humanoid Robot Kangaroo

Humanoid bipedal platforms are finally moving from controlled lab environments to real-world applications in unstructured scenarios. This has been recently achieved with substantial improvements in hardware and software architectures. In particular, the hardware design of bipedal humanoid robots has changed to be robust and resilient to impacts, as well as lightweight and efficient for dynamic movements. At the same time, planning and control architectures permit now to generate and control high dynamic, contact rich, motions employing the dynamics and the full body of the system.

In this short paper, we present Kangaroo, a new humanoid bipedal robot designed by PAL Robotics for research on agile and dynamic locomotion. Kangaroo exploits novel mechanical solutions based on linear actuators and closed/parallel kinematics chains. We present the mechanical design and the whole planning and control pipeline developed to make Kangaroo jump on real experiments to demonstrate the superior robustness and resilience of the platform.

I. INTRODUCTION

Boston Dynamics videos had demonstrate that humanoid robots could outperform human capabilities, by realizing fast and dynamic motions that most humans are not able to execute. To achieve those skills is necessary a perfect combination of planning and robust control software able to react against disturbances, a proper hardware design that offers enough speed and torque to perform those motions, and a reliable platform with high resilience to impacts. Following this objective, PAL Robotics has designed the new bipedal research platform Kangaroo for exploring advanced control methods for legged locomotion 1 .

The robot, with 6 actuated Degrees of Freedom (DOFs) per leg, has been designed to be lightweight (12.5 kg per leg) with low moving inertia so that it is suitable for highly dynamic motions such as jumping and running. The design of Kangaroo takes advantage of closed linkage mechanisms to place all the actuators close to the torso, protecting the motors and electronics against impacts and shaping the gear ratio curve to provide high speed or high torque depending on the configuration.

Kangaroo also introduces a novel kinematic design in which the swing and the leg extension/contraction motions {adria.roig, sai.kishor, narcis.miguel, enrico.mingo, luca.marchionni}@pal-robotics.com ‡ Toward, Boulevard Lazare Carnot 81, 31000 Toulouse, France. E-mail: {pierre.fernbach}@toward.fr 1 https://pal-robotics.com/robots/kangaroo/ are completely decoupled. To jump vertically or to perform vertical motions, one single actuator per leg is activated, being highly efficient compared to many other robots in which this movement is achieved by the simultaneous activation of multiple actuators at the hip, knee, and ankle, acting in opposite directions.

All the actuators are driven with a brushless motor connected to a ball screw, providing high efficiency, low friction, and resistance to high forces in the linear direction. In the hip and ankle roll and pitch joints, linear actuators are arranged in differential configurations to increase the limb's force capability without substantially increasing reflected inertia. In addition, Kangaroo employs a new generation of in-house electronics boards for motor control and sensor acquisition, that allows embedding custom and decentralized closed-loop force/torque controllers at the actuator level. This permits running low-level controllers at a higher frequency, up to 20 kHz. The actuators are equipped with force sensors that perceive the external forces applied through the mechanism.

In terms of software, Kangaroo uses ROS as a standard framework [START_REF] Quigley | Ros: an open-source robot operating system[END_REF] benefiting from the available functionalities for testing, visualization, and communication between processes, making it a versatile platform for researchers. ros control [START_REF] Chitta | ros control: A generic and simple control framework for ros[END_REF] is in charge to deal with hardware and software abstraction layers, thereby exposing the interfaces to the controllers. This short paper introduces the platform, the main software architecture, and preliminary results achieved with Kangaroo.

II. RELATED WORKS

Several bipedal humanoid robots present closed linkages, in particular parallel four-bars mechanisms for the actuation of the pitch-roll ankle joints [START_REF] Knabe | Team valor's escher: A novel electromechanical biped for the darpa robotics challenge[END_REF]- [START_REF] Radford | Valkyrie: Nasa's first bipedal humanoid robot[END_REF]. However, in this case, the actuators are moved along the shin, making the advantages of such designs localized only to the ankle part.

The Atrias robot [START_REF] Hubicki | Atrias: Design and validation of a tether-free 3d-capable spring-mass bipedal robot[END_REF] is the first of a series of famous biped systems which uses closed linkages on the legs kinematics. The new BostonDynamics Atlas [START_REF] Guizzo | By leaps and bounds: An exclusive look at how boston dynamics is redefining robot agility[END_REF] and AgilityRobotics Cassie/Digit [START_REF] Apgar | Fast online trajectory optimization for the bipedal robot cassie[END_REF], [START_REF] Hurst | Walk this way: To be useful around people, robots need to learn how to move like we do[END_REF], both presents closed linkages in their mechanical design. In particular, Cassie/Digit uses electric actuators localized near its torso/pelvis with a transmission, consisting of closed-bar linkages, to transfer the motion to reduce inertia, especially during impacts while walking. Finally, RH5v2 is a recently developed, full humanoid robot, presenting closed kinematic loops both in the upper and lower body [START_REF] Mronga | Whole-body control of seriesparallel hybrid robots[END_REF].

Following this new trend in the design of humanoid robots with closed and parallel kinematic linkages, this work presents the preliminary results on the modeling, planning, and control of the lower body of the Kangaroo platform.

III. THE KANGAROO ROBOT

The Kangaroo robot has been designed using novel highpower and robust linear actuation units located near the pelvis area, while the motion is transferred to the joints through a complex system of closed and parallel linkages. Kangaroo is actuated by 12 linear ball-screw modules with electrical motors, 6 per leg, while the complete leg can be modeled by 76 DOFs, among them, 64 are passive. Figure 1 shows the location of the actuators in Kangaroo, 5 out of 6 are placed at the hip and only 1 is placed at the femur. Table I reports a comparison between linear actuators used in robotics and the 2 types of units used in Kangaroo 2 .

The actuation is reported at the joint-level by a series of closed and parallel mechanisms consisting of several four-bar linkages, differential linkages, and other closed-bar linkages.

The peculiar design of these linkages on Kangaroo permits to extend the feet by maintaining constant the orientation of the ankle w.r.t. the base frame emulating a prismatic joint, and also to move the hip in the roll, pitch, and yaw direction, as well to move the ankle in the roll and pitch direction. In particular, the extension/retraction of the leg is achieved by just actuating the motor on the femur. Notice that the hip pitch and roll, as well the ankle pitch and roll are coupled, by the respective parallel linkages, forming complex nonlinear differential joints.

We developed a reduced model to curtail the complexity of the full model with 76 DOFs and to avoid the use of closed kinematic chains that are not fully supported by many simulators or kinematics/dynamics libraries. This reduced model of Kangaroo is used to perform Operational Space Inverse model. We employ a virtual linear joint between the hip and the ankle to emulate the effect of the transmission for the leg and the second open kinematic branch contains the femur and the tibia links starting at the hip and terminating at the ankle. Since this formulation contains two open chains, they need to be constrained inside the inverse dynamics (see Section V). Finally, the quantities computed using this model are then mapped to actuation space using transmissions which considers the full system of closed linkages (see Section IV).

IV. ACTUATION MAPPING THROUGH TRANSMISSION

As previously written, Kangaroo control software infrastructure is based on ros control, which provides a set of powerful libraries and tools that helps to define the structure and a properly defined functionality of each of its components. One such component in this framework is the Transmissions.

Transmissions interfaces are the point where the joint commands are transformed into the actuator commands, and vice versa. ros control provides already the interface for the most common transmissions such as revolute joints, prismatic joints, or differential joints.

At this point, the transmission component of the ROS framework helps us to build the Kangaroo custom nonlinear transmission plugins based on the geometrical system of the leg where a single actuator moves multiple passive joints. This makes sure that the complexity of the closed-loop linkages is pushed down to the low-level components and thereby exposing a simplified joint structure control system, the one of the reduced model, to the high-level components such as the OSID controller.

Kangaroo has 4 different custom nonlinear transmission plugins built to control the whole leg structure:

• Virtual leg length joint, • Ankle pitch and roll joint, • Hip roll and pitch joint,

• Hip yaw joint. One of the benefits to hide all the closed linkage complexity in the transmission component is the re-usability of previously developed controllers since most of them command the robot in joint space, and not in actuation space.

V. QP-BASED OPERATIONAL SPACE INVERSE DYNAMICS

The operational space tasks are transformed into torque references using an OSID stack solved using Quadratic Programming (QP) optimization. In addition to the classical formulation of the OSID problem, we consider as well the presence of a closed linkage in our model. Being q ∈ R 16 + SO(3) and ν ∈ R 22 the robot generalized coordinates and velocities, respectively, the equations of motions can be expressed as:

M(q) ν + h(q, ν) = Sτ + J T c (q)F c + J T v (q)λ, (1)
with F c ∈ R 24 and λ ∈ R 4 , constrained by

J c (q) ν + Jc (q, ν)ν = 0, (2)
accounting for the contacts and

J v (q) ν + Jv (q, ν)ν = 0, (3)
accounting for the closed linkages between the hip and the ankles. Notice that S considers both the underactuation associated with the floating base and the open linkages in the model. We designed a novel OSID scheme based on three main steps:

• First, a QP constrained by the first 6 rows of ( 1) is solved to find optimal ν * and F * c , considering for instance, tasks such as the Center of Mass (COM), angular momentum, base orientation and other constraints such as friction cones.

• Secondly, the Lagrange multipliers λ associated to the closed linkages are solved exploiting the particular structure of the reduced model3 :

λ * = J -T v,u (q) (M j,u (q) ν * + h j,u (q, ν)) , (4)
where the subscripts u and j mean the underactuated and joint space part, respectively, of the matrix where the subscripts are applied. • Third, the actuated torques τ are retrieved substituting at the actuated part of (1) the obtained ν * , F * c and λ * .

This OSID algorithm is integrated as a ROS controller that sends computed optimal torques to a decentralized, high frequency, low-level torque controller, before sending the commands to the actuators. A further component consists of a kinematic state estimator for the pose and the velocities of the base of the robot, needed by the OSID. Finally, the references for the OSID are computed as a plugin within the controller.

This architecture allows sharing the same structure among different behaviors: stabilizing, walking, jumping, etc.

VI. JUMP PLANNING AND STABILIZATION

The jumping controller generates the references to make the robot jump, rectify the position and orientation in the air, and stabilize it after landing. Since the objectives and constraints involved in each phase of the jump are different they need to be replaced online.

During the first phase of the jumping, the desired height of the jump is transformed into a ramp velocity profile of the COM with a constant acceleration that will depend on the desired height. In this first stage, the classical contact constraints, i.e. friction cones and feet contacts, are active. During this phase, it is important to minimize the angular momentum to avoid undesired rotations during the next flight phase.

Once the COM reaches a desired velocity or height, the robot is considered to be in the flight phase. The flight phase consists of a parabolic shot in which the time of flight and COM height will depend on the last velocity achieved at the end of the first phase. At this point, all the contact constraints have to be deactivated and only the postural task and the angular momentum task will remain active. The angular momentum task helps to stabilize and correct possible rotations due to non-zero angular velocities at the moment of jumping.

When contacts are detected, the landing phase starts. In this phase, we add a constraint to control desired resultant ground reaction forces for each foot and the friction cone constraints and the COM tasks are activated again. The reference ground reaction forces for the landing phase are computed off-line employing trajectory-optimization based on CasADi [START_REF] Andersson | CasADi -A software framework for nonlinear optimization and optimal control[END_REF]. The forces go from high, to resist the impact, to low, to actively absorb the landing, back to the nominal value to keep the height of the COM (see Figure 3). During the first and final phases, the whole motion is stabilized through a task at the

ẍ = ω 2 (x -z d ) , (5)
where x and ẍ are the projected position and accelerations of the COM, respectively, z d is the reference position of the ZMP, and ω = g h . z d is computed instantaneously from the desired capture point (CP) positions and velocities, in fact:

c = x + 1 ω ẋ, (6a)
→ ∂ ∂t ċ = ẋ + 1 ω ẍ = ẋ + ω(x -z), → z d = c - 1 ω ċd , (6b)
with c and ċ the CP and its velocity, respectively. The term ċd can be treated as a new reference input based on a P feedback law plus feed-forward:

ċd = ċr + K(c r -c). (7)

VII. RESULTS

We tested our preliminary control architecture on the Kangaroo robot to perform a series of jumps. The robot successfully jumped at various heights several times demonstrating robustness and repeatability, even though it was not able to always land correctly without falling, especially for higher jumps where the landed COM position stays far from the center of the support polygon. Fails in landings were mostly caused by issues related to force mapping from joints to actuators and state-estimation. Figure 4 shows one of the 10 cm jumps with a correct landing, on the top sequence, and a higher jump, more than 20 cm, with a failed landing, on the bottom sequence. The jumping videos can be found at https://www.youtube.com/ playlist?list=PLg5jtQ hM1RezuPS6CCQpbiv27L ZAdtI.

VIII. CONCLUSION AND FUTURE WORKS

In this short paper, we presented the Kangaroo platform, introducing its mechanical design, and the planning and control software pipeline. The entry point of this pipeline is the reduced model which simplifies the complex kinematic In the close future, we plan is to substitute this, so-called jumping controller, with a model predictive controller, for example, based on the centroidal dynamics of the system. This will improve disturbances rejection, particularly for higher jumps. Future works also include to feed-forward the trajectories generated through the existing trajectory optimization framework to generate different kinds of motions such as lateral jumps, forward jumps, spins, etc. Finally, the employment of the full model for the OSID computation will also be considered to accurately compensate for the dynamics of the system.

Fig. 1 :

 1 Fig. 1: Side and front view of the Kangaroo robot and actuators placing.

2Fig. 2 :

 2 Fig. 2: On the left, the kinematics of the reduced model with the open and closed branch. On the right the reduced model visualized in RVIZ, again with open and closed branch.

Fig. 3 :

 3 Fig.3: Contact forces at the right foot for the jump, planned using trajectory optimization.

  a. Successful landing of 10 cm jump b. Failure landing of more than 20 cm jump

Fig. 4 :

 4 Fig. 4: In the upper sequence Kangaroo jumping 10 cm and successfully landing. In the lower sequence a high jump with a failed landing.

TABLE I :

 I Comparison between linear actuators

		Apptronik	PAL	Moog-IIT ISA [12]	RH5 Knee [13]
		P170 Orion [11]
			LA5	LA2	v2	v5
	Peak Force [kN ] Max Speed [ m s ] Mass [kg]	3.2 0.27 0.786	5.0 0.65 2.1	2.0 0.42 0.78	4.0 -0.920	7.5 -1.6	5.845 0.14 -
	Stroke [cm]	6.47	15.0	8.0	8.0	10.0	11.8

In this particular case, the contact forces F * c does not appear in the computation


Sarah Temmam

Khamsing Vongphayloth

Eduard Baquero

Sandie Munier

Massimiliano Bonomi

Béatrice Regnault

Bounsavane Douangboubpha

Yasaman Karami

Delphine Chrétien

Daosavanh Sanamxay

Eduard Baquero Salazar

Vilakhan Xayaphet

Phetphoumin Paphaphanh

Vincent Lacoste

Somphavanh Somlor

Khaithong Lakeomany

Nothasin Phommavanh

Philippe Pérot

Océane Dehan

Faustine Amara

Flora Donati

Thomas Bigot

Michael Nilges

Félix A Rey

Sylvie Van Der Werf

Paul T Brey

Marc Eloit
email: marc.eloit@pasteur.fr

Bat coronaviruses related to SARS-CoV-2 and infectious for human cells

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

Introduction

The origin of SARS-CoV-2, as its mode of introduction into the human population, is currently unknown. Since its emergence, numerous animal species have been studied to identify possible reservoirs and/or intermediate hosts of the virus, including a large diversity of insectivorous bats of the genus Rhinolophus. Despite the recent report of various SARS-CoV-2-related viruses in Rhinolophus shameli (isolated in Cambodia in 2010 [START_REF] Delaune | A novel SARS-CoV-2 related coronavirus in bats from Cambodia[END_REF] ), R. pusillus and R. malayanus (China, 2020 and 2019 respectively 2 ), in R. acuminatus (Thailand, 2020 3 ) and R. cornutus (Japan, 2013 4 ), the closest SARS-CoV-2 bat-borne genome still remains the one from R. affinis, RaTG13 (China, 2013) [START_REF] Zhou | A pneumonia outbreak associated with a new coronavirus of probable bat origin[END_REF][START_REF] Rahalkar | Lethal Pneumonia Cases in Mojiang Miners (2012) and the Mineshaft Could Provide Important Clues to the Origin of SARS-CoV-2[END_REF] with 96.1% similarity at the whole genome level. Several studies also suggested the involvement of pangolin coronaviruses in the emergence of SARS-CoV-2 [START_REF] Liu | Are pangolins the intermediate host of the 2019 novel coronavirus (SARS-CoV-2)?[END_REF][START_REF] Xiao | Isolation of SARS-CoV-2-related coronavirus from Malayan pangolins[END_REF][START_REF] Wahba | An Extensive Meta-Metagenomic Search Identifies SARS-CoV-2-Homologous Sequences in Pangolin Lung Viromes[END_REF] . Since its appearance in humans, SARS-CoV-2 has evolved through sporadic mutations and recombination events [START_REF] Jackson | Generation and transmission of interlineage recombinants in the SARS-CoV-2 pandemic[END_REF] , some of which correspond to gains in fitness allowing the virus to spread more widely, or to escape neutralizing antibodies [START_REF] Rochman | Ongoing global and regional adaptive evolution of SARS-CoV-2[END_REF] .

To decipher the origin of SARS-CoV-2, it is therefore essential to ascertain the diversity of animal coronaviruses and more specifically of bat coronaviruses. Although the identification of SARS-CoV-2 in bats is a major goal, a more realistic objective is to identify the sequences that contribute to its mosaicism. The spike sequence appears essential, as it determines the binding affinity and accessibility of the receptor-binding domain (RBD) to the cellular angiotensin-converting enzyme 2 (ACE2) receptor and is therefore responsible for host range [START_REF] Letko | Functional assessment of cell entry and receptor usage for SARS-CoV-2 and other lineage B betacoronaviruses[END_REF][START_REF] Shang | Structural basis of receptor recognition by SARS-CoV-2[END_REF][START_REF] Wang | Structural and Functional Basis of SARS-CoV-2 Entry by Using Human ACE2[END_REF] . The closest related bat strain identified to date (RaTG13) has a low RBD sequence similarity to SARS-CoV-2 and, with only 11/17 human ACE2 (hACE2) contact amino-acid residues conserved with SARS-CoV-2, its affinity to hACE2 is very limited [START_REF] Liu | Binding and molecular basis of the bat coronavirus RaTG13 virus to ACE2 in humans and other species[END_REF] . Moreover, SARS-CoV-2 poorly infects bats and bat cells tested so far [START_REF] Aicher | Species-specific molecular barriers to SARS-CoV-2 replication in bat cells[END_REF] . In addition, no bat SARS-CoV-2-like virus has been shown to use hACE2 to efficiently enter human cells, and none presents the furin cleavage site that is associated with an increased pathogenicity in humans [START_REF] Johnson | Loss of furin cleavage site attenuates SARS-CoV-2 pathogenesis[END_REF] . SARS-CoV-2 RBD binds to Rhinolophus macrotis ACE2 with a lower affinity than to hACE2 [START_REF] Liu | Cross-species recognition of SARS-CoV-2 to bat ACE2[END_REF] . An essential piece of information, i.e. finding bat viruses with an RBD motif close to that of SARS-CoV-2 and capable of binding to hACE2 with high affinity, is therefore missing.

We hypothesized that this type of virus could be identified in bats living in the limestone karstic terrain common to China, Laos, and Vietnam in the Indochinese peninsula. We report here the presence of sarbecoviruses close to SARS-CoV-2 whose RBDs differ from that of SARS-CoV-2 by only one or two contact residues, that strongly bind to the hACE2 protein, and mediate a hACE2-dependent entry and replication into human cells. Despite the absence of the furin cleavage site, these viruses may have contributed to SARS-CoV-2's origin and may intrinsically pose a future risk of direct transmission to humans.

Main

Diversity of bat and coronavirus species

A total of 645 bats belonging to six families and 46 species were captured (Table S1). Two hundred and forty-seven blood samples, 608 saliva, 539 anal/feces, and 157 urine swabs were collected from the northern part of Laos (Table S2). We first screened all 539 feces samples through a pan-coronavirus RT nested PCR approach [START_REF] Chu | Avian Coronavirus in Wild Aquatic Birds[END_REF] . Overall, 24 individuals of 10 species were positive, and one individual (BANAL-27) was concomitantly infected by an alphacoronavirus and a betacoronavirus (Table S3). BLAST analysis of amplicons identified alphacoronavirus sequences of the Decacovirus, Pedacovirus, and Rhinacovirus subgenera and betacoronavirus sequences of the Nobecovirus and Sarbecovirus subgenera. Sequences of the Sarbecovirus subgenus were all identified from Rhinolophus individuals belonging to three different species i.e., R. malayanus, R. marshalli, and R. pusillus. Positive individuals were trapped in three different districts, and those infected with a sarbecovirus were all from the Fueng district in Vientiane Province (Fig. 1A, site 1).

NGS together with Sanger sequencing was used to obtain a complete genomic sequence of five of the seven sarbecoviruses (Fig. 1, Table S4). The coverage of the genome of the remaining two sarbecoviruses, i.e. BANAL-27 and BANAL-242 sampled from R. pusillus and R. malayanus bats respectively, was 90% and therefore they were not included in the final analyses. Phylogenetic analyses performed on the receptor-binding protein domain of lineages A and B human SARS-CoV-2 [START_REF] Rambaut | A dynamic nomenclature proposal for SARS-CoV-2 lineages to assist genomic epidemiology[END_REF] , and on representative bat and pangolin sarbecoviruses, placed the Laotian R. malayanus BANAL-52, R. pusillus BANAL-103, and R. marshalli BANAL-236 coronaviruses close to human SARS-CoV-2 and pangolin coronaviruses collected in 2019, while R. malayanus BANAL-116 and BANAL-247 coronaviruses belonged to a sister clade with other bat coronaviruses (RmYN02, RacCS203, RpYN06, and PrC31) from different Rhinolophus species. Pangolin coronaviruses sampled in 2017 displayed a basal position relative to these strains (Fig. 1B). Interestingly, one should note that very similar SARS-CoV-2-like viruses are shared by different bat species, suggesting a possible circulation of viruses between different species living sympatrically in the same caves. These results are consistent with the similarity plot analysis showing that RaTG13 and BANAL-52 bat coronaviruses exhibit high nucleotide identity with SARS-CoV-2 throughout the length of the genome (96.8% for BANAL-52 and 96.1% for RaTG13). Interestingly, BANAL-52 presents a higher level of nucleotide conservation than RaTG13 in the S1 domain of the spike, and especially in the spike's N-terminal domain (NTD) and RBD (Fig. 1C). These observations are congruent with amino-acid identities between human SARS-CoV-2 and representative bat and pangolin coronaviruses, which present a high level of conservation, except for ORF8 of bat BANAL-116, BANAL-247, Rc-o319 and RmYN02 (Supp. Figure 1). Interestingly, the S1 domain of the spike (and especially the N-terminal domain) presents a lower degree of conservation in several bat coronaviruses, suggesting that this domain may reflect a relative degree of adaptation of the virus to its mammalian host (Fig. 1D, Extended data 1).

Bat Sarbecovirus evolutionary history

Following GARD analysis, we identified 14 recombinant breakpoints during the evolutionary history of sarbecoviruses, which were further confirmed by phylogenetic analyses performed on the 15 fragments of sequences defined by the breakpoints (Fig. 2, Supp. Figure 2). No specific signature was identified in the breakpoints (Table S5). SARS-CoV-2 presents a mosaic genome, to which contributed more than five sequences close to sequences published or determined during this study: R. malayanus RmYN02 and R. pusillus RpYN06 viruses found in China in 2019, R. affinis RaTG13 coronavirus found in China in 2013, and R. malayanus BANAL-52 and R. pusillus BANAL-103 found in North Laos in 2020 (this study). No pangolin coronavirus sequence was immediately associated with a recombination event at the origin of SARS-CoV-2. Laotian Rhinolophus bat coronaviruses presented a lower degree of recombination compared to SARS-CoV-2. Such recombination events occurred between other BANAL viruses isolated from bats living sympatrically in caves in the same area.

Interestingly, the origin of several fragments of SARS-CoV-2 genomes could be assigned to several donor strains and not a unique donor sequence. For example, a breakpoint was identified 7 nucleotides upstream of the RBD region of S1: the downstream fragment of SARS-CoV-2, which comprises the RBD and the beginning of S2, could involve BANAL-52 R. malayanus, BANAL-103 R. pusillus, and BANAL-236 R. marshalli viruses, which formed a highly supported sister clade of SARS-CoV-2 (fragment 11, Supp. Figure 2). In a more basal position are R. shameli bat coronaviruses and pangolin-2019 coronaviruses, consistent with the conservation of RBD amino-acid sequences among SARS-CoV-2 and representative bat and pangolin coronaviruses (Extended data 2). Among the 17 residues that interact with human ACE2, 16 are conserved between SARS-CoV-2 and BANAL-52 or -103 (one mismatch, H498Q), and 15/17 are conserved for BANAL-236 (two mismatches, K493Q and H498Q) while only 13/17 residues are conserved for the Cambodian bat R. shameli virus and 11/17 for the Chinese bat R. affinis RaTG13 virus. At the full spike protein level, bat R. affinis RaTG13 and pangolin-2017 P4L viruses looked closer to SARS-CoV-2 than bat R. malayanus BANAL-52, but this effect is due to a higher degree of conservation in S2. All these viruses shared the absence of a furin cleavage site and the conservation of the internal fusion peptide (Supp. Figures 1 &3 and Extended data 3).

Interaction of BANAL RBDs with ACE2

Biolayer interferometry experiments to measure the interaction between hACE2 and the RBDs of BANAL-52/103 (which have identical residues in the receptor-binding motif) (Extended data 2), BANAL-236 and SARS-CoV-2 (residues 233-524), resulted in a dissociation constant Kd three times lower for the BANAL RBDs compared to SARS-CoV-2 (Fig. 3A). This higher affinity can be attributed to the Q498H mismatch that has been reported to increase the affinity of SARS-CoV-2 RBD for hACE2, and also to be involved in the host range expansion of SARS-CoV-2 and SARS-CoV-2 like viruses [START_REF] Wrapp | Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation[END_REF][START_REF] Laffeber | Experimental Evidence for Enhanced Receptor Binding by Rapidly Spreading SARS-CoV-2 Variants[END_REF][START_REF] Lei | Neutralization of SARS-CoV-2 spike pseudotyped virus by recombinant ACE2-Ig[END_REF][START_REF] Walls | Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein[END_REF][START_REF] Niu | Molecular basis of cross-species ACE2 interactions with SARS-CoV-2-like viruses of pangolin origin[END_REF][START_REF] Zhang | SARS-CoV-2 rapidly adapts in aged BALB/c mice and induces typical pneumonia[END_REF][START_REF] Huang | Q493K and Q498H substitutions in Spike promote adaptation of SARS-CoV-2 in mice[END_REF][START_REF] Zhang | Bat and pangolin coronavirus spike glycoprotein structures provide insights into SARS-CoV-2 evolution[END_REF] .

To study the effect of the mutations at the interface between these RBDs and hACE2, we performed Molecular Dynamics (MD) simulations of the SARS-CoV-2, BANAL-236 and BANAL-52/103 RBD/hACE2 complexes initiated from the crystal structure and homology models of these systems, respectively (Table S6). Cluster analysis of the MD trajectories revealed that, at the RBD-hACE2 interface, both BANAL complexes were identical to the SARS-CoV-2 RBD/hACE2 complex within 2 Å backbone RMSD (Extended data 4), except for one of the BANAL-52/103 simulations that displayed larger conformational variability of the RBD residues S443-Y449 (Supp. Methods and Extended data 5). Empirical scoring functions predicted a similar RBD/hACE2 binding energy in all three complexes (Extended data 6).

The analysis of the persistence of hydrogen bonds and salt bridges provided further insights into the effect of the substitutions at the RBD/hACE2 interface (Fig. 3B and Extended data 7). The H498Q mismatch present in both BANAL-52/103 and BANAL-236 RBDs disrupted the hydrogen bonds between RBD Q498 and both hACE2 K353 and Q42. However, these hydrogen bonds were only transiently formed in the SARS-CoV-2 complex: more persistent hydrogen bonds in this region (RBD T500 -hACE2 D355, RBD G502 -hACE2 K353, RBD Y505 -hACE2 E37) were not affected. The K493Q mismatch enabled the formation of two salt bridges between RBD and hACE2 that were not present in the SARS-CoV-2 complex (RBD K493 -hACE2 E35 and RBD K493 -hACE2 D38).

For further insight into the molecular details of these interactions, we determined the crystal structure of the complex BANAL-236 RBD with the hACE2 peptidase domain to 2.9 Å resolution (Table S7). The overall structure of this RBD is identical to that of the SARS-CoV-2 (RMSD 0.360 Å, 150 C). The only significant difference is in the region between amino acids D363 and S375 (Fig. 3C,arrow). In this region, BANAL spikes have the A372T mismatch, which converts the sequence 370 NSA 372 to 370 NST 372 , an N-linked glycosylation sequon. The crystals indeed showed clear electron density for the first Nacetylglucosamine residue of glycan attached to N 370 . The main chain at residue T 372 makes a hydrogen bond with the glycan residue, altering the conformation of the main chain downstream which results in partial unwinding of helix H4 in RBD 236. The calculation of a simulated annealing composite omit map for the segment D363-S375 confirms the correct assignment of the structure for this polypeptide segment and the N-glycosylation at residue N370 (Extended data 8).

As expected, most of the interactions observed in the SARS-CoV-2 RBD/hACE2 complex [START_REF] Lan | Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor[END_REF] are also present in the structure of the BANAL-236 RBD/hACE2 complex. In these interfaces, there are three main clusters of interactions as indicated in Fig. 3C (insets). The sequence mismatches are in clusters 2 (making the salt bridge RBD K493 -hACE2 E35) and 3 (hydrogen bond between RBD H498 and D38). Although the interaction K493-E35 contributes to stabilizing the complex, it does not seem to affect drastically the binding to hACE2 because both BANAL RBDs have similar Kd values.

Virus replication in human cells

To assess whether the BANAL-236 spike protein could mediate entry into cells expressing human ACE2, we generated lentiviral particles pseudotyped with the Wuhan or the BANAL-236 spikes (Supp. Figure 4). We detected spike-mediated entry of the BANAL-236 spike-pseudotyped lentivirus in 293T-ACE2, contrarily to control cells not expressing hACE2 (Fig. 4A). Entry was blocked by human sera neutralizing SARS-CoV-2, but not by control non-neutralizing sera, demonstrating that neutralization of BANAL-236 was specific for epitopes shared with the spike of SARS-CoV-2 (Fig. 4B).

In order to isolate infectious viruses, rectal swabs were inoculated on VeroE6 cells. No CPE was observed 3 and 4 days after infection, but viral RNAs were detected for one of the 2 wells inoculated with the BANAL-236 sample (CT = 25.1 at D3, CT = 21.7 at D4). The culture supernatant (C1) formed plaques on VeroE6 and the titer was 3800 pfu/mL. A C2 viral stock was prepared by amplification on VeroE6 at a MOI of 10 -4 . The culture supernatant was harvested on day 4 when CPE was observed and titrated on VeroE6 (Extended data 9). The plaques' phenotype was small, but the titer reached 2.6.10 6 pfu/mL. The random NGS performed on the RNA extracted from this stock confirmed that the culture was pure and corresponded to the BANAL-236 virus, without any non-synonymous mutations between the original BANAL-236 genome and the C2 viral stock. Replication of BANAL-236 in VeroE6 was efficiently inhibited by soluble human ACE2, thus showing that entry and propagation was largely ACE2-dependent (Fig. 4C, Extended data 9). Furthermore, BANAL-236 replicated in human cell lines expressing endogenous levels of ACE2, Calu-3 and Caco-2 (Fig. 4C). The kinetics of RNA synthesis were slower compared to Wuhan virus. Infectious viral particles were produced at D4 in the supernatant of Caco-2 and Calu-3 cells (respectively 10 4.7 and 10 2.9 pfu/mL).

Discussion

Many sarbecoviruses circulate in Rhinolophus colonies living in caves in China and probably also in neighboring countries further south [START_REF] Hu | Discovery of a rich gene pool of bat SARS-related coronaviruses provides new insights into the origin of SARS coronavirus[END_REF][START_REF] Ge | Isolation and characterization of a bat SARS-like coronavirus that uses the ACE2 receptor[END_REF][START_REF] Latinne | Origin and cross-species transmission of bat coronaviruses in China[END_REF] . During the course of a prospective study in Northern Laos, we have identified, among other coronaviruses, five sarbecoviruses for which we obtained full-length sequences. Among these, three (BANAL-52, -103, and -236) were considered to be close to SARS-CoV-2 because of the similarity of one of the S1 domains (NTD, RBD) or S2 with SARS-CoV-2.

Because genomic regions subject to recombination are likely contributing to host-virus interactions, we compared SARS-CoV-2 strains from the two lineages identified at the very onset of the COVID-19 outbreak [START_REF] Rambaut | A dynamic nomenclature proposal for SARS-CoV-2 lineages to assist genomic epidemiology[END_REF] to these novel bat sarbecoviruses and to pangolin strains within the SARS-CoV-2 clade. Strains close to R. pusillus RpYN06, R. malayanus RmYN02, and Rhinolophus sp. PrC31 isolated in China in 2018-2019, along with R. malayanus BANAL-52, R. pusillus BANAL-103, and R. marshalli BANAL-236 isolated in Laos in 2020, contributed to the appearance of SARS-CoV-2 in different regions of the genome. No closer viral genome has yet been identified as a possible contributor, and pangolin coronaviruses appear as more distantly related than bat coronaviruses. We identified potential recombination sites, allowing for the reconstruction of the phylogenetic history of early isolated SARS-CoV-2 strains between homologous regions defined by recombination points. We identified a breakpoint at the beginning of the SARS-CoV-2 RBD, resulting in a downstream fragment key for the virus tropism and host spectrum composed of the RBD, the furin cleavage site, and ending in the Nterminal region of S2. Despite the absence of the furin site, phylogenetic reconstruction of this fragment revealed that Laotian R. malayanus BANAL-52, R. pusillus BANAL-103, and R. marshalli BANAL-236 coronaviruses are the closest ancestors of SARS-CoV-2 known to date. ORF8 was highly divergent between SARS-CoV-2 related genomes. ORF8 from strains BANAL-52, -103, -236, like that of RaTG13, were closer to SARS-CoV-2 than to pangolin strains. ORF8 encodes a protein that has been proposed to participate in immune evasion [START_REF]Novel Immunoglobulin Domain Proteins Provide Insights into Evolution and Pathogenesis of SARS-CoV-2-Related Viruses[END_REF] and is deleted in many human SARS-CoV-2 strains that appeared after March 2020 [START_REF] Su | Discovery and Genomic Characterization of a 382-Nucleotide Deletion in ORF7b and ORF8 during the Early Evolution of SARS-CoV-2[END_REF] , which is reminiscent of the deletions identified during the 2003 SARS epidemic [START_REF] Chinese | Molecular evolution of the SARS coronavirus during the course of the SARS epidemic in China[END_REF] . Therefore, the presence of ORF8 is consistent with bats acting as a natural reservoir of early strains of SARS-CoV-2.

Structural and functional biology studies have identified the RBD domain that mediates the interaction with hACE2 and host range, as well as the major amino acids that are involved [START_REF] Lan | Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor[END_REF][START_REF] Conceicao | The SARS-CoV-2 Spike protein has a broad tropism for mammalian ACE2 proteins[END_REF][START_REF] Damas | Broad host range of SARS-CoV-2 predicted by comparative and structural analysis of ACE2 in vertebrates[END_REF] . The RBDs (BANAL-52, -103, and -236) are closer to SARS-CoV-2 than that of any other bat strain described so far, in particular that of RaTG13, the virus identified in R. affinis from the Mojiang mineshaft where pneumonia cases with clinical characteristics a posteriori interpreted as similar to COVID-19 [START_REF] Rahalkar | Lethal Pneumonia Cases in Mojiang Miners (2012) and the Mineshaft Could Provide Important Clues to the Origin of SARS-CoV-2[END_REF] were recorded in 2012 [START_REF] Cohen | Wuhan coronavirus hunter Shi Zhengli speaks out[END_REF][START_REF] Ge | Coexistence of multiple coronaviruses in several bat colonies in an abandoned mineshaft[END_REF] . Overall, one (H498Q (BANAL-103 and -52)) or two (K493Q and H498Q (BANAL-236)) amino acids interacting with hACE2 are substituted in these strains in comparison to SARS-CoV-2. These mutations did not destabilize the BANAL-236 / hACE2 interface, as shown by the BLI experiments (Fig. 3A) and analyzed by MD simulations.

Our results therefore support the hypothesis that SARS-CoV-2 could originally result from a recombination of sequences pre-existing in Rhinolophus bats living in the extensive limestone cave systems of South-East Asia and South China [START_REF] Clements | Limestone Karsts of Southeast Asia: Imperiled Arks of Biodiversity[END_REF][START_REF] Hassanin | Inferring the ecological niche of bat viruses closely related to SARS-CoV-2 using phylogeographic analyses of Rhinolophus species[END_REF] . Many species forage in the same cave areas, including R. malayanus and R. pusillus [START_REF] Soisook | A taxonomic review of Rhinolophus stheno and R. malayanus (Chiroptera: Rhinolophidae) from continental Southeast Asia: an evaluation of echolocation call frequency in discriminating between cryptic species[END_REF] . In addition, the distribution of R. marshalli, R. malayanus, and R. pusillus overlaps in the Indochinese sub-region (Supp. Figure 5), which means they may share caves as roost sites and foraging habitats [START_REF] Francis | Field Guide to the Mammals of South-east Asia[END_REF] . With the novel viruses here described, understanding the emergence of SARS-CoV-2 does not need to hypothesize recombination or natural selection for increased RBD affinity for hACE2 in an intermediate host like the pangolin before spillover [START_REF] Makarenkov | Horizontal gene transfer and recombination analysis of SARS-CoV-2 genes helps discover its close relatives and shed light on its origin[END_REF] , nor natural selection in humans following spillover [START_REF] Andersen | The proximal origin of SARS-CoV-2[END_REF] . However, we found no furin cleavage site in any of these viruses on sequences determined from original fecal swab samples, devoid of any risk of counterselection of the furin site by amplification in Vero cells [START_REF] Johnson | Loss of furin cleavage site attenuates SARS-CoV-2 pathogenesis[END_REF] . Lack of furin cleavage site may be explained by insufficient sampling in bats. Based on comparison of the sequences around the S1/S2 cleavage site (Extended data 3), it has been suggested that the furin cleavage site present in SARS-CoV-2 could originate from recombination events between SARS-CoV-2-related coronaviruses co-circulating in bats 2,[START_REF]The Sarbecovirus origin of SARS-CoV-2's furin cleavage site -SARS-CoV-2 coronavirus / nCoV-2019 Evolutionary History[END_REF] , meaning that BANAL-116, BANAL-247, bat RmYN02 2 and bat RacCS203 3 coronaviruses may share a common history with SARS-CoV-2. Alternatively, the furin cleavage site could have been acquired through passages of the virus in an alternate host or during an early poorly symptomatic unreported circulation in humans. Finally, the epidemiological link between these bat viruses and the first human cases remains to be established.

As expected from the high affinity for ACE2 of the S ectodomain of BANAL-236, pseudoviruses expressing it were able to enter efficiently human cells expressing endogenous hACE2 using an ACE2dependent pathway. However, alternative routes of entry may still exist, especially in cells that do not express ACE2 [START_REF] Puray-Chavez | Systematic analysis of SARS-CoV-2 infection of an ACE2-negative human airway cell[END_REF] . Entry was blocked by a serum neutralizing SARS-CoV-2. The RaTG13 strain, the closest to SARS-CoV-2 known before, had never been isolated. In contrast, preliminary studies show that BANAL-236 replicated in primate VeroE6 cells with a small plaque phenotype compared to SARS-CoV-2. Further analysis may indicate more clearly which steps shape infectivity.

To conclude, our results pinpoint the presence of new bat sarbecoviruses that seem to have the same potential for infecting humans as early strains of SARS-CoV-2. Guano collectors, or certain ascetic religious communities who spend time in or very close to caves, as well as tourists visiting caves, are particularly at risk of being exposed. Further investigations are needed to assess if such exposed populations have been infected, symptomatically or not, by one of these viruses, and whether infection could confer protection against subsequent SARS-CoV-2 infections. In this context, it is noteworthy that SARS-CoV-2 with the furin site deleted replicates in hamsters and in transgenic mice expressing hACE2, but leads to less severe disease while protecting from rechallenge with wild-type SARS-CoV-2 [START_REF] Johnson | Loss of furin cleavage site attenuates SARS-CoV-2 pathogenesis[END_REF] .

Figure 1. Genomic description of bat-borne sarbecoviruses identified in Laos. (A)

Map of sampling sites. All BANAL isolates were collected from the same site (site 1). Map was downloaded from https://www.divagis.org/gdata. (B) Phylogenetic analysis of the protein sequence of the receptor-binding domain of Laotian and representative human, bat, and pangolin sarbecoviruses. Sequences were aligned with MAFFT [START_REF] Katoh | MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization[END_REF] in "auto" mode, and maximum likelihood phylogenetic reconstruction was performed with PhyML implemented through the NGPhylogeny portal [START_REF] Lemoine | fr: new generation phylogenetic services for non-specialists[END_REF] with the LG+G substitution model. Branch support was evaluated with the aBayes parameter. Accession numbers and bat species are specified in the name of the sequences. Sequences are colored according to Fig. 1C. (C) Similarity plot analysis of Laotian and representative bat and pangolin sarbecoviruses based on the full-length genome sequence of SARS-CoV-2 human prototype strain (NC_045512, Wuhan-Hu-1) used as reference. The analysis was performed with the Kimura-2 parameter model, a window size of 1,000 base pairs, and a step size of 100 base pairs with SimPlot program, version 3.5.1 [START_REF] Lole | Full-Length Human Immunodeficiency Virus Type 1 Genomes from Subtype C-Infected Seroconverters in India, with Evidence of Intersubtype Recombination[END_REF] . (D) Heatmap of identities at the protein level of representative human, bat, and pangolin sarbecoviruses compared to human SARS-CoV-2 lineage B (NC_045512). Spike protein has been divided into functional domains, and the sequences are ordered according to percentage of identity of the RBD domain. "*": absence of a functional ORF10 in Thai bat RacCS203 (accession number MW251308). Heatmap was created using the gplots package in R (version 3.6.3). Extended data 1. Spike identity matrices at the genus level of representative sarbecoviruses. Amino-acid (lower) and nucleotide (upper) identity matrices of Laotian and representative human, bat, and pangolin sarbecoviruses. Spike N-terminal (NTD), Receptor-binding (RBD) and S2 nucleotide and amino-acid sequences were aligned with MAFFT, and identity matrices were constructed using CLC Main Workbench 21.0.4 (Qiagen). Matrices were colored according to the identity scale, from 25% (red) to 100% (green) of nucleotide or aminoacid identity. Extended data 4. Analysis of the stability and conformational heterogeneity of RBD/hACE2 complexes. Time series (left column) and violin plots (right column) of backbone Root Mean Square Deviation (RMSD) from the initial, energy-minimized model calculated on the residues in RBD (A), hACE2 (B), at the interface of RBD and hACE2 (C), and on the entire complex (D). In the violin plots, the white circle corresponds to the median value, the black rectangle extends from the first to the third quantiles, and the thin black line represents the 95% confidence intervals. Population of the 3 most significant clusters visited during the course of the MD simulations (E). The analysis is performed for 9 different MD simulations: 3 replicates of the SARS-CoV-2 (shades of green), BANAL-236 (shades of red), and BANAL-52/103 (shades of blue) RBD/hACE2 complexes. Extended data 6. Estimation of RBD-hACE2 binding energy. Time series (left column) and violin plots (right column) of the RBD-hACE2 binding energy estimated using ROSETTA (A) and FoldX (B). In the violin plots, the white circle corresponds to the median value, the black rectangle extends from the first to the third quantiles, and the thin black line represents the 95% confidence intervals. The analysis is performed for 9 different MD simulations: 3 replicates of the SARS-CoV-2 (shades of green), BANAL-236 (shades of red), and BANAL-52/103 (shades of blue) RBD/hACE2 complexes.

Extended

Extended data 5. Additional analysis of the MD

Extended data 7. Analysis of the inter-subunits hydrogen bonds at

Methods

Ethical and legal statements

The bat study was approved by the wildlife authorities of the Department of Forest Resource Management (DFRM), and the Ministry of Agriculture and Forestry, Lao PDR, No. 2493/DFRM, issued on May 21, 2020; No. 0755/MAF issued on June 2, 2020. All animals were captured, handled, and sampled following previously published protocols and ASM guidelines [START_REF]Protocol for Bat and Rodent Sampling Methods[END_REF][START_REF] Sikes | Guidelines of the American Society of Mammalogists for the use of wild mammals in research[END_REF] . Exportation from Lao and importation in France were conducted according to national regulations. Human serum samples used for neutralization assays were already available [START_REF] Virachith | Low seroprevalence of COVID-19 in Lao PDR, late 2020[END_REF] and selected based on their status regarding seroneutralization of SARS-CoV-2. They were collected following a protocol approved by the Lao National Ethics Committee for Health Research (NECHR) (Ref #052/2020).

Biosafety

Both the Institut Pasteur du Laos (IPL) and the Faculty of Environmental Science have extensive experience in safely collecting bats (appropriate biosafety training and Personal Protective Equipment (PPE) for collectors). In this study, training prior to field work was organized once for field work participants. The aim was to teach participants on the transmission risks of infectious agents from bats and how to identify, assess, and mitigate these risks, as well to practice the use of PPE.

During field collection, sampling stations were selected to minimize potential exposure to infectious agents and stress on the animals during the handling time by selecting (1) an area easy to disinfect, (2) out of view of the general public, (3) a location that would not expose the general population (like a picnic area), and (4) procedures that reduce time and stress on bats caused by handling. For handling bats and bat samples the following minimum PPE was required: eye protection, an N95 respirator, long clothing/cover all, and latex gloves (2pairs). All wastes were disposed of in biohazard bags and were transported to a disposal site of IPL in Vientiane capital. Each sample box containing 81 cryovial tubes was put inside individual plastic ziplock and transported from field to IPL using dry ice by packing them in a polystyrene foam box. At IPL, samples were stored in a specific -80°C freezer until analysis.

For initial sample analysis at IPL, samples were transferred to BSL-3 room for nucleic acid extraction. Initial Coronavirus screening by RT-PCR was then performed under internal regulation on biosafety and security of IPL. Sample extracted products (50 µl) were stored in NucleoSpin 8 sample boxes (8x12) and sent to IP-Paris by packing them in the individual plastic ziplock and using dry ice. Aliquots of anal swab samples were sent in a separated box by individually triple packaging in dry ice.

All experiments on potentially infectious samples done at Institut Pasteur (Paris) were conducted in BSL-3 laboratories according to procedures adapted to respiratory viruses.

Bat sampling areas and sample collection

Trapping sessions were conducted on four sites, in Fueng and Meth Districts, Vientiane Province, and in Namor and Xay Districts, Oudomxay Province, between July 2020 and January 2021 (Fig. 1A, Tables S1 andS2). Bats were captured using four-bank harp traps [START_REF] Francis | A Comparison of Mist Nets and Two Designs of Harp Traps for Capturing Bats[END_REF] and mist nets set in forest patches between rice fields/orange/banana plantations and karst limestone formations, for 5-8 nights depending on accessibility. Harp traps were set across natural trails in patches of forest understory. Mist nets were set across natural trails, at the edges of forests, at entrances of caves, and in areas near cave entrances, as well as in open areas or those with high forest canopy. Bats were morphologically identified following morphological criteria [START_REF] Francis | A Comparison of Mist Nets and Two Designs of Harp Traps for Capturing Bats[END_REF][START_REF] Hutson | Mammals of the Indomalayan Region: A Systematic Review by[END_REF][START_REF] Csorba | Horseshoe Bats of the World (Chiroptera : rhinolophidae). Alana Books Bishop's Castle Shropsh[END_REF] . Other data such as forearm length (FA), sex, developmental stage (adult or juvenile) and reproductive condition (pregnant or lactating) were also recorded. Bats were sampled for saliva, feces and/or urine, and blood before release at the capture site. Species identification of PCR-positive individuals was confirmed by sequencing the mitochondrial Cytochrome oxidase 1 [START_REF] Ge | Isolation and characterization of a bat SARS-like coronavirus that uses the ACE2 receptor[END_REF] .

Initial Coronavirus screening

Total RNA was extracted from feces samples using the NucleoSpin 8 virus kit (Macherey Nagel). cDNA was synthesized using the Maxima H minus first strand cDNA synthesis kit (Thermo Scientific) and random hexamers following the manufacturer's instructions. The presence of coronaviruses was tested by a nested PCR approach using PCR master mix (Promega) and by targeting the RNA-dependent RNA Polymerase gene using combinations of degenerate and non-degenerate consensus primers as previously described [START_REF] Chu | Avian Coronavirus in Wild Aquatic Birds[END_REF] . PCR products of the expected size were directly sequenced on both strands by Sanger sequencing using the nested PCR primers. The sequences obtained were confirmed by similarity analysis using the NCBI BLASTn search (http://www.ncbi.nlm.nih.gov/BLAST).

Primer design for Betacoronavirus enrichment before Next-Generation Sequencing

Betacoronavirus enrichment was performed at the genus level by adapting a previously described protocol [START_REF] Deng | Metagenomic sequencing with spiked primer enrichment for viral diagnostics and genomic surveillance[END_REF] based on k-mers for targeted-sequence enrichment prior to NGS. Briefly, 2,000 complete Betacoronavirus genomes were downloaded from the GenBank and GISAID databases and then clustered to a 95% sequence identity using CD-HIT-EST [START_REF] Fu | CD-HIT: accelerated for clustering the next-generation sequencing data[END_REF] . Overall, 185 representative sequences of all betacoronaviruses were used for further analysis. Due to the high diversity within the genus Betacoronavirus, the full genomes belonging to the subgenera Sarbecovirus, Nobecovirus, and Merbecovirus plus Embecovirus were separately aligned using MAFFT multiple sequence alignment software and used to design 13-mer spiked primers per cluster. The genomic position of the full set of primers was extracted from some representative subgenera sequences and close primers were removed. Finally, 416 spiked primers were synthesized by Eurofins Genomics Germany GmbH.

Sample preparation for sequencing

Reverse transcription was performed using the mix of spiked primers and random hexamers at a 10:1 ratio using the SuperScript™ IV First-Strand Synthesis System (Invitrogen). After a denaturation step at 95°C for 3 min in the presence of dNTPs (500 µM), DTT (5 mM), and RNaseOUT inhibitor, the RT reaction was incubated for 10 thermal cycles consisting of 6 steps at 8°C for 12 s, 15°C for 45 s, 20°C for 45 s, 30°C for 30 s, 35°C for 2 min and 42°C for 3 min, followed by a final incubation step at 42°C for 20 min, as previously described [START_REF] Regnault | Deep Impact of Random Amplification and Library Construction Methods on Viral Metagenomics Results[END_REF] . Double-stranded cDNA was generated using the Sequenase2.0 DNA Polymerase (Applied Biosystems™) in presence of dNTPs and then purified using the Beckman Coulter™ Agencourt AMPure XP.

For samples with a low nucleic acid content, a random amplification step was performed using the MALBAC Single Cell WGA kit (Yikon Genomics, Promega). The amplified product was then purified using AMPure XP Beads, eluted in a final volume of 20 μL of low TE (10 mM Tris-HCl (pH 8,0), 0,1 mM EDTA), and quantified with the Qubit™ DNA HS Assay (Life Technologies, Thermo Fisher Scientific Inc.).

Libraries were generated using the NEBNext Ultra II DNA Library Prep kit (New England Biolabs) after a fragmentation step using the Covaris M220 Focused-ultrasonicator using microTUBE-15 (Peak Incident power (W)=18, Duty Factor=20%, Cycles per Burst=50, Treatment time (sec)=60). The PCR amplified libraries were cleaned up using 0.9X AMPure XP Beads and checked on the 2100 Bioanalyzer system with the High-Sensitivity DNA kit (Agilent Technologies, France) and quantified with the Qubit™ DNA HS Assay. Finally, the dual-multiplexed libraries were pooled (six samples/pool) and run on the Illumina NextSeq500 platform with High Output Kit v2.5 (150 Cycles).

Amplicon sequencing

In addition to enrichment-based sequencing, cDNA was amplified using the AmpliSeq for Illumina SARS-CoV-2 Research Panel (cat# 20020496), applying twenty-six amplification cycles in PCR1 and nine cycles in PCR2. Primers at the end of the amplicons were partially digested during the library preparation, following manufacturer's instructions. Libraries were barcoded individually using the Illumina UD dual-indexes and normalized with the AmpliSeq Library Equalizer for Illumina (cat# 20019171), then pooled and sequenced on the Illumina NextSeq500 instrument using a Mid Output v2.5 kit (SR 150 cycles). To eliminate residual PCR primer sequences, raw reads were trimmed by 15 bases at each end and special attention was paid to checking that internal sequences corresponding to primer regions in overlapping amplicons did not derive from primer sequences of the multiplex PCR. As an internal control, we verified that sequences of the complete genome of sample BANAL-236 obtained from the enrichment-based sequencing approach and from the AmpliSeq approach were identical.

Genome assembly and finishing

Raw reads from the enrichment-based sequencing were processed with an in-house bioinformatics pipeline (Microseek; Bigot et al., submitted) comprising quality check and trimming, reads normalization, de novo assembly, and ORF prediction of contigs and singletons, followed by 3 levels of taxonomic assignation [START_REF] Bratuleanu | The virome of Rhipicephalus, Dermacentor and Haemaphysalis ticks from Eastern Romania includes novel viruses with potential relevance for public health[END_REF] . Sequences identified as Sarbecovirus were then mapped onto appropriate reference sequences using CLC Genomics Workbench 20.0 (Qiagen). Trimmed reads from the amplicon sequencing were mapped to the SARS-CoV-2 genome first, then mapped again (refined mapping) to the closest genome relative. When needed, complete genomes were obtained by conventional PCR and Sanger sequencing. Briefly, viral RNA was reverse transcribed using the SuperScript IV reverse transcriptase (Invitrogen, USA) and cDNA was subsequently used to fill the gaps in the genomes using the Phusion High Fidelity DNA polymerase (New England Biolabs, France) and specific primers flanking the missing regions. Positive PCR products were further purified and sequenced by Sanger sequencing at Eurofins Genomics.

Recombination and phylogenetic analyses

Identification of recombination events occurring during the evolutionary history of bat sarbecoviruses was performed using the IDPlot package [START_REF] Pond | Automated Phylogenetic Detection of Recombination Using a Genetic Algorithm[END_REF] , a web-based workflow that includes multiple sequence alignment and phylogeny-based breakpoints prediction using the GARD algorithm from the HyPhy genetic analysis suite [START_REF] Pond | GARD: a genetic algorithm for recombination detection[END_REF] . First, a comprehensive analysis comprising 106 sequences and covering all nonhuman Sarbecovirus and Sarbecovirus-related complete genomes available in GenBank and GISAID databases was performed, including prototype strains of SARS-CoV-2 isolated in 2019. Special attention was paid to including bat-borne and pangolin-borne viral sequences in order to maximize the ability to capture a large diversity of the sarbecoviruses. Then, a reduced set of 36 sequences was chosen because of their phylogenetic proximity with SARS-CoV-2, and the GARD algorithm was run to identify recombination breakpoints in Laotian and representative human, bat, and pangolin sarbecoviruses. Breakpoint coordinates were confirmed by performing phylogenetic analyses on the corresponding fragments using PhyML implemented through the NGPhylogeny portal [START_REF] Lemoine | fr: new generation phylogenetic services for non-specialists[END_REF] . Branch support was evaluated with the aBayes parameter.

Generation of lentiviral pseudoviruses

The BANAL-236 and Wuhan synthetic spike genes were cloned into the pVAX1 vector with a cytoplasmic tail truncation of 19 amino acids. Pseudotyped lentiviral particles were prepared using HEK-293T cells (ATCC CRL-3216) seeded in 10-cm dishes. 293T cells were co-transfected with a 5 µg of spike-encoding plasmid, 10 µg of lentiviral backbone plasmid expressing the firefly luciferase (pHAGE-CMV-Luc2-IRES-ZsGreen-W), and 3.3 µg of each lentiviral helper plasmid expressing HIV Gag-Pol (HDM-Hgpm2), Tat (HDM-tat1b), and Rev (pRC-CMV-Rev1b) using calcium-phosphate precipitation [START_REF] Crawford | Protocol and Reagents for Pseudotyping Lentiviral Particles with SARS-CoV-2 Spike Protein for Neutralization Assays[END_REF] . The medium was replaced 5 h post-transfection by 6 mL of DMEM without fetal calf serum (FCS) and phenol red. Pseudotyped particles were harvested 48 h post-transfection, clarified by centrifugation at 2500 g for 5 min and frozen at -80°C. Mock pseudotyped lentivirus was generated as above but in the absence of an S-expressing plasmid.

Spike-pseudotyped lentivirus entry assays

HEK-293T cells stably expressing human ACE2 (293T-ACE2) were transduced in suspension by mixing 50 µL of 3-fold serial dilutions of S-pseudotyped lentiviruses with 50 µL of cells at 4.10 5 cells/mL in 96well white culture plates [START_REF] Anna | High seroprevalence but short-lived immune response to SARS-CoV-2 infection in Paris[END_REF] . At 60-72 h post-transduction, 100 µL of Bright-Glo luciferase substrate (Promega) were added to the wells and luminescence was measured using a Berthold Centro XS luminometer.

Neutralization assays

Sera neutralizing or not SARS-CoV-2 were described in [START_REF] Virachith | Low seroprevalence of COVID-19 in Lao PDR, late 2020[END_REF] . They were decomplemented at 56°C for 30 min and 2.5 µL were incubated with 0.5 µL of S-pseudotyped lentiviruses in a final volume of 50 µL of DMEM-10% FCS without phenol red in 96-well white culture plates. After 30 min at room temperature, 50 µL of 293T-ACE2 cells in suspension at 4.10 5 cells/mL were mixed to the wells. Luminescence was measured at 60-72 h post-transduction as described above. Neutralization was calculated using the following formula: 1 -(RLU in presence of serum / (mean of RLU in absence of serum determined in 12 wells -3 x STD).

Virus isolation and multiplication

Rectal swabs were inoculated in duplicate in 24-well plates containing VeroE6 cells (ATCC CRL-1586) (1/5 dilution in 100 µL of DMEM without FCS supplemented with 1% Penicillin-Streptomycin, 1% Fungizone, and 1 µg/mL TPCK-treated trypsin). After 1 h of adsorption at 37°C, inoculum was removed and 1 mL of the medium described above was added. Three and 4 days after infection, cytopathic effect (CPE) was monitored and 100 µL of supernatant was collected for RNA extraction. RT-qPCR targeting a conserved sequence in the E gene was performed as described [START_REF] Corman | Detection of 2019 novel coronavirus (2019-nCoV) by real-time RT-PCR[END_REF] . Culture supernatant (C1) was harvested at day 4 and titrated by plaque assay on VeroE6 overlaid with 0.5% carboxymethylcellulose containing 1 µg/mL TPCK-treated trypsin. A viral stock was prepared by amplification on VeroE6 cells at a MOI of 10 -4 . Culture supernatant (C2) was harvested at day 4 when massive CPE was observed and titrated by plaque assay on VeroE6 as described above. A C3 viral stock was produced by a subsequent viral amplification at a MOI of 10 -4 for 3 days on VeroE6. RNA was extracted from the viral stock and submitted to random NGS analysis using the SMARTer® Stranded Total RNA-Seq Kit v3 -Pico Input Mammalian (Takara Bio USA, Inc.). Raw reads were processed with Microseek pipeline, as described above.

For infection experiments, VeroE6 or human cell lines Calu-3 (lung cells, ATCC HTB-55) and Caco-2 (intestinal cells, ATCC HTB-37) were infected in triplicate in 24-well plates with BANAL-236 (C3) or Wuhan viral stocks at a MOI of 0.0001 and 0.01 respectively. Human soluble ACE2 was pre-incubated at 25 µg/mL 67 for 30 min with the viral inoculum before VeroE6 infection. Infections were carried out without TPCK in the medium described above for 3 and 4 days. Supernatants were recovered at 0-, 3and 4-days post-infection and genome copy number and viral titers were quantified as described above.

Protein expression and purification

BANAL-52/103, BANAL-236 and SARS-CoV-2 RBDs (residues 233-524, with C-terminal 8xHis-Strep and AVI tags) and hACE2 peptidase domain (residues 19-615, with C-terminal 8xHis tag) were expressed in Expi293F cells at 37°C and 8% CO2 (GnTI-Expi293™, ThermoFisher Scientific). Cell culture supernatants were collected five days post transfection and purified by affinity chromatography followed by size exclusion chromatography (SEC) using a 200 10/300 GL column pre-equilibrated in 20 mM Tris-HCl pH 8.0, 100 mM NaCl.

For crystallization experiments, the same constructs were expressed in Expi293F GnTI cells. The protein tags were cleaved overnight with thrombin and deglycosylated with EndoH. The RBD was mixed with a 1.3 molar excess of hACE2 and the complex was purified by SEC.

Biolayer Interferometry

Purified Avi-tagged RBD was biotinylated using a BirA biotin-protein ligase kit according to manufacturer's instructions (Avidity). The biotinylated RBDs at 100 nM were immobilized to SA sensors. A 1:2 dilution series of hACE2 starting at 100 nM in PBS-BSA buffer was used in cycles of 200 s association followed by 200 s dissociation steps to determine protein-protein affinity. The data were baseline-subtracted and the plots fitted using the Pall FortéBio/Sartorius analysis software (version 12.0). Data were plotted in Prism 9.1.0.

Crystallization and data collection

Crystals of complex BANAL-236 RBD/hACE2 were obtained at 4°C in sitting drops by mixing 200 nl of the protein complex at 8 mg/mL with 200 nl of reservoir solution containing 0.2 M lithium sulfate, 0.1 M Tris 8.5, 30 % w/v PEG 4000. The crystals were soaked in reservoir solution containing 20% glycerol as cryoprotectant before being flash-frozen in liquid nitrogen. X-ray diffraction data was collected on the beamline PROXIMA 1 at the SOLEIL synchrotron (St Aubin, France) and reduced using the XDS package 68 . The structure of the complex was determined by molecular replacement with Phaser software [START_REF] Mccoy | Phaser crystallographic software[END_REF] using the coordinates of SARS-CoV-2 RBD in complex with hACE2 as search template (Protein Data Bank (PDB) 6M0J). The model was manually corrected in COOT [START_REF] Emsley | Coot: model-building tools for molecular graphics[END_REF] and refined with phenix.refine [START_REF] Liebschner | Macromolecular structure determination using X-rays, neutrons and electrons: recent developments in Phenix[END_REF] . The final coordinates were deposited in the PDB with the entry code 7PKI.

Molecular Dynamics simulations of RBD/hACE2 complexes

Generation of homology models. Homology models of the BANAL-236 and BANAL-52/103 RBD/hACE2 complexes were constructed using the X-ray structure of the SARS-CoV-2 RBD/hACE2 complex (PDB code 6M0J; resolution 2.45 Å) as template using MODELLER v. 10.1 [START_REF] Sali | Comparative protein modelling by satisfaction of spatial restraints[END_REF] . These models included the fragments S19-D615 and T333-G526 of hACE2 and RBD, respectively, which were the regions resolved in the template. In these regions, BANAL-236 and BANAL-52/103 RBDs have a sequence identity to SARS-CoV-2 RBD equal to 96.9% and 97.4%, respectively. The alignment reported in Extended data 3 was used. The zinc and chloride atoms present in the template were retained during homology modeling; N-Acetylglucosamine (NAG) and water residues were removed. 7 disulfide bonds were detected by MODELLER in the template (3 in hACE2 and 4 in the RBD) and enforced in the generation of the homology models using CHARMM-like distance and dihedral angles restraints. For each construct, 100 homology models were built and ranked based on the normalized DOPE score [START_REF] Shen | Statistical potential for assessment and prediction of protein structures[END_REF] . The top 3 scoring models of each complex were used as starting points of 3 independent MD simulations, as described in the following section.

Details of the MD simulations: setup, equilibration, and production. The 6 homology models described in the previous section along with the X-ray structure of the SARS-CoV-2 RBD/hACE2 complex (PDB code 6M0J) were used as input to the CHARMM-GUI server [START_REF] Jo | CHARMM-GUI: a web-based graphical user interface for CHARMM[END_REF] . In the case of the X-ray structure, for consistency with the homology models, the zinc and chloride atoms present in 6M0J were retained while the NAG and water residues were removed. The 7 systems were solvated in a triclinic box of initial x-y-z dimensions of ~13.5 nm * 9.2 nm * 8.3 nm. Potassium and chloride ions were added to ensure charge neutrality at a salt concentration of 0.15 M. The total number of atoms was ~104000. Additional details of the systems are reported in Table S4. The CHARMM36m force field [START_REF] Huang | CHARMM36m: an improved force field for folded and intrinsically disordered proteins[END_REF] was used for the protein and ions and the TIP3P model [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF] was used for the water molecules. CHARMM36m force field parameters for the 7 pairs of cysteines linked by disulfide bonds were used. The CHARMM-GUI models were first energy-minimized using the steepest descent algorithm. After minimization, the systems were equilibrated using a 1-ns long simulation in the NPT ensemble followed by a 1-ns long simulation in the NVT ensemble. The temperature T was set at 300 K and the pressure P at 1 atm using the Bussi-Donadio-Parrinello thermostat [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] and the Berendsen barostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] , respectively. During equilibration, harmonic restraints on the positions of the protein backbone and sidechains heavy atoms were applied. For each system studied, production simulations were performed in the NVT ensemble for 1 µs. A time step of 2 fs was used together with LINCS constraints on h-bonds [START_REF] Hess | LINCS: A linear constraint solver for molecular simulations[END_REF] . The van der Waals interactions were gradually switched off at 1.0 nm and cut off at 1.2 nm; the particlemesh Ewald method was used to calculate electrostatic interactions with cutoff at 1.2 nm [START_REF] Essmann | A smooth particle mesh Ewald method[END_REF] . Production simulations were performed at room temperature for consistency with the BLI experiments.

Details of the analysis. To evaluate the stability of the starting model during the production simulations, we calculated the backbone Root Mean Square Deviation (RMSD) with respect to the energy-minimized structure for each frame of the trajectories. The RMSD was calculated separately for the residues in: RBD (Extended data 5A), hACE2 (Extended data 5B), RBD-hACE2 interface (Extended data 5C), and in the entire complex (Extended data 5D). Interfacial residues were defined as the residues in one subunit closer than 0.8 nm to the residues in the other subunit in the X-ray structure of the SARS-Cov-2 RBD/hACE2 complex. RMSD calculations were performed using the driver utility of PLUMED v. 2.7 [START_REF]Promoting transparency and reproducibility in enhanced molecular simulations[END_REF] . To characterize conformational heterogeneity at the RBD-hACE2 interface, we performed a cluster analysis of the MD trajectories using as similarity metrics the backbone RMSD of the interfacial residues and the GROMOS clustering approach [START_REF] Daura | Peptide folding: When simulation meets experiment[END_REF] with a cutoff equal to 0.2 nm. The 9 trajectories were first concatenated, clustering was then performed, and finally the population of each cluster was calculated separately for each trajectory (Extended data 5E). To estimate the binding energy between RBD and hACE2, we used the InterfaceAnalyzer tool in ROSETTA v. 3.11 [START_REF] Stranges | A comparison of successful and failed protein interface designs highlights the challenges of designing buried hydrogen bonds[END_REF] (Extended data 6A) and the AnalyseComplex tool in FoldX v. 4 [START_REF] Schymkowitz | The FoldX web server: an online force field[END_REF] (Extended data 6B). To identify relevant interactions at the RBD and hACE2 interface that could contribute to the binding affinity of the complex, we quantified the frequency of formation of inter-subunits salt bridges (Fig. 3B) and hydrogen bonds Extended data 7) during the course of the MD simulations. For each frame of the trajectories, we used PLUMED to calculate the distances between the sidechain charged groups of aspartic acids (OD1/OD2), glutamic acids (OE1/OE2), lysines (NZ), and arginines (NH1/NH2). An inter-subunits salt bridge was then defined as formed if the distance between groups with opposite charge was lower than 0.32 nm. We confirmed this calculation using the Salt Bridges tool available in VMD [START_REF] Humphrey | visual molecular dynamics[END_REF] . To monitor the formation of inter-subunits hydrogen bonds, we used the Hydrogen Bond Analysis module of the MDAnalysis library v. 1.0.0 [START_REF] Michaud-Agrawal | MDAnalysis: a toolkit for the analysis of molecular dynamics simulations[END_REF] . A donor-acceptor distance and angular cutoffs of 0.3 nm and 150° were used to define the formation of a hydrogen bond. We confirmed this calculation using the Hydrogen Bonds tool available in VMD. Data processing was carried out using XDS [START_REF] Monteil | Inhibition of SARS-CoV-2 Infections in Engineered Human Tissues Using Clinical-Grade Soluble Human ACE2[END_REF] . Molecular replacement and model refinement were performed using PHENIX 68,[START_REF] Emsley | Coot: model-building tools for molecular graphics[END_REF] . Model building was carried out using COOT [START_REF] Mccoy | Phaser crystallographic software[END_REF] . Structural figures were made using PYMOL (The PyMOL Molecular Graphics System, v.2.0 Schrödinger). We used 2 crystals for X-ray data collection. Each crystal resulted in one set of X-ray data. The best dataset (as judged by data statistics) was used for structure determination and refinement.

*Values in parentheses are for the highest-resolution shell.

Table S8: Acknowledgment of sharing of SARS-CoV-2-related genome sequences from the GISAID database related to Figures 1-2 and Supp. Figures 1234. We thank the authors listed in the table below for sharing their genomic sequences of coronaviruses analyzed in this study. The identification of the breakpoints is highly dependent on the programs used, the sequences datasets and the quality of the alignments, leading to difficulties in comparing independent studies performed with different tools. Phylogenetic analyses performed on different fragments are therefore necessary to confirm the identification of breakpoints, especially by confirming that the position of a given sequence differs between fragments bracketing the breakpoint. This is the strategy presented below: sequences are colored according to Figure 2.

Accession

The substitution model considered for all phylogenetic reconstructions presented below was GTR+G. PhyML was used to infer phylogenies both at the whole genome and at the breakpoints nucleotide level. The tree topology search was performed by selecting the best of NNI (Nearest Neighbor Interchange) and SPR (Subtree Pruning and Regraphing). The statistical test selected for branch support was the approximate Bayes parameter, as proposed in the NGPhylogeny portal. Only posterior probabilities >0.5 were considered.

Figure 2 .

 2 Figure 2. Recombination events in the evolutionary history of sarbecoviruses. Schematic representation of the 15 recombinant fragments of relevant Sarbecovirus genomes compared to SARS-CoV-2 human prototype strain (NC_045512). Breakpoint positions are referred to the position in the alignment. Where possible, the closest viral sequence is indicated for each fragment. In other cases, "MULT" (group of multiple sequences) is mentioned. "$": unresolved fragment phylogeny (fragment 13, from positions 27344 to 27800 in the alignment). Sequences are colored as in Fig. 1. The complete phylogenetic analyses are presented in Supp.Figure 2.

  Figure 2. Recombination events in the evolutionary history of sarbecoviruses. Schematic representation of the 15 recombinant fragments of relevant Sarbecovirus genomes compared to SARS-CoV-2 human prototype strain (NC_045512). Breakpoint positions are referred to the position in the alignment. Where possible, the closest viral sequence is indicated for each fragment. In other cases, "MULT" (group of multiple sequences) is mentioned. "$": unresolved fragment phylogeny (fragment 13, from positions 27344 to 27800 in the alignment). Sequences are colored as in Fig. 1. The complete phylogenetic analyses are presented in Supp.Figure 2.

Figure 3 .

 3 Figure 3. Dynamics of the binding of hACE2 to bat-sarbecovirus borne RBDs and structural insight of the complex. (A) Biolayer interferometry binding analysis of the hACE2 peptidase domain to immobilized BANAL52/103, BANAL-236 RBDs or SARS-CoV-2. Black lines correspond to global fit of the data using a 1:1 binding model. (B) Frequency of formation of salt bridges close to the RBD/ACE2 interface (from left to right: D30/K417, E35/K493, D38/K493, K31/E35, and D38/K353) during the course of the MD simulations. The analysis is performed for 9 different MD simulations (3 replicates for each complex) of hACE2 in complex with SARS-CoV-2 (shades of green), BANAL-236 (shades of red) and BANAL-52/103 RBDs (shades of blue). (C) Ribbon representations of the crystal structures of hACE2 peptidase domain (cyan) in complex with SARS-CoV-2 (PDB 6M0J) or BANAL-236 (this study, PDB 7PKI) RBDs (pink). Black arrows in the overall structures indicate the structural difference between the two complexes at the level of helix H4. The insets show the main interactions in the ACE2-RBD interfaces. Residues in the RBM mutated between SARS-CoV-2 and BANAL-236 are indicated with colored boxes.

Figure 4 .

 4 Figure 4. BANAL-236 entry and propagation in human cells. (A) Results of spike-pseudotyped BANAL-236 (square) and Wuhan-Hu-1 (diamond) pseudoviruses entry assay in HEK-293T cells expressing (purple) or not (grey) the hACE2 receptor, expressed in Relative Luminescence unit (RLU) produced by the firefly luciferase present in the lentiviral backbone and the Bright-Glo luciferase substrate. A single experiment performed in triplicate representative of two experiments is shown. Center values represent the average of the 3 replicates and error bars indicate SD. (B) Results of spike-pseudotyped BANAL-236 (black) and Wuhan-Hu-1 (grey) neutralization assay expressed in % of neutralization of luciferase activity in the absence of serum. SARS-CoV-2 neutralizing sera were from patients with confirmed infections while non-neutralizing sera samples were collected before SARS-CoV-2 spread in Laos. Dashed line: neutralization threshold. A single experiment representative of three independent experiments is shown. (C) Human cell lines expressing endogenous ACE2, Calu-3 (blue line) and Caco-2 (green line), were infected at a MOI of 0.01 with BANAL-236 (square dots) and Wuhan (diamond dots) viruses. VeroE6 cells (red line) were infected at a MOI of 0.0001 with BANAL-236 and Wuhan viruses pre-incubated or not (grey line) with soluble human ACE2 (sACE2) at 25 µg/mL for 30 min. Genome copy number was quantified by RT-qPCR in the supernatants recovered 3-and 4-days post-infection. A single experiment performed in triplicate is shown. Center values represent the average of the 3 replicates and error bars indicate SD.

data 2 .

 2 Alignment of spike RBD domain. Protein alignment of the Receptor Binding Domain (RBD) of Laotian and representative human, bat and pangolin sarbecoviruses. Sequences were aligned with MAFFT in G-iNS-I mode. Residues interacting with human ACE2 receptor are highlighted in grey. The domain used for interactions modeling, based on the X-ray structure 6M0J (residues T333 to G526), is highlighted by a black line. Extended data 3. Nucleotide and amino-acid alignments of the furin cleavage site region. Complete nucleotide and amino-acid spike sequences of representative bat SARS-CoV-2-like coronaviruses were downloaded from GenBank and GISAID and aligned with MAFFT (G-INS-I parameter) (A & C). Alignments were manually edited as proposed by Zhou 2 and Lytras 47 with CLC Main Workbench (Qiagen) (B & D). Alignments of the furin cleavage region are presented at the nucleotide (A & B) and the amino-acid (C & D) level, respectively.

  simulations of the BANAL-52/103 RBD/hACE2 complex. Comparison of the time series of interface RMSD during the course of two MD simulations of the BANAL-52/103 RBD/hACE2 complex with short (BANAL-52/103-CoV.1) and long (BANAL-52/103-CoV.1*) equilibration phase. The large fluctuations of the interface RMSD are due to the flexibility of the RBD loop between residues S443 and Y449 (insets, in yellow). When these residues were not included in the calculation of the interface RMSD, the time series displayed a more stable behavior (BANAL-52/103-CoV.1-L and BANAL-52/103-CoV.1*-L).

  the interface of RBD and hACE2. Frequency of formation of hydrogen bonds at the interface of RBD and hACE2 in the knob (A), base (B), and tip regions (C). The analysis is performed for 9 different MD simulations: 3 replicates of the SARS-CoV-2 (shades of green), BANAL-236 (shades of red), and BANAL-52/103 (shades of blue) RBD/hACE2 complexes. Extended data 8. Sticks representation of segments D364-S375 of BANAL 236 (left panel) and SARS-CoV-2 (right panel) RBDs. A 2Fo-Fc composite omit map (contoured at 3is shown for this region in BANAL 236 RBD. Extended data 9. Isolation of BANAL-236 on VeroE6 cells. (A) CPE observed on VeroE6 4 days after inoculation at an MOI of 10-4 from the C1. (B) Uninfected VeroE6 cells layer. (C) Plaque assay performed from the C2 stock on VeroE6 cells. (D) Comparative CPE observed on VeroE6 cells infected with Wuhan (top panel) or BANAL-236 (bottom panel) in absence (left) or in presence (right) of soluble ACE2. Photos were taken with a EVOS XL Core microscope at x10. A single experiment performed in triplicate is shown for each cell line.

Figure S2 :

 S2 Figure S2: Complete phylogenetic reconstructions of recombinant fragments, as described in Figure 2.

Figure S3 :

 S3 Figure S3: Alignment of the whole spike amino-acid sequences.

Figure S4 :

 S4 Figure S4: Western-Blot analysis of concentrated lentivral particles pseudotyped with Wuhan or BANAL-236 spike proteins.Pseudotyped lentivral particles were concentrated by ultracentrifugation at 40,000g for 2h through a 20% sucrose cushion (MLS-50 rotor in an Optima MAX-XP ultracentrifuge). Pellets were resuspended in 250 µL of DMEM for 30 min at 4°C. An equal amount of 2X NuPAGE LDS sample buffer was added and different volumes were loaded onto a 4-12% NuPAGE Bis-Tris gel (Thermofisher). PVDF membranes (Hybond®, Amersham) were incubated with primary antibodies directed against Spike (ABIN199984, 1/2,000) or p24 (NIH, 1/2000) and revealed with peroxidase-conjugated secondary antibodies (GE Healthcare) and the ECL2 substrate (Thermofisher). The chemiluminescence signals were acquired using a ChemiDoc Imager (Bio-Rad) (uncropped blot below). A single experiment was performed, and cropped and un-cropped blots are presented below.

Figure S5 :

 S5 Figure S5: Distribution of bats species from which SARS-CoV-2-like viruses were detected in this study.Distribution of R. pusillus (green), R. marshalli (red), and R. malayanus (blue). Data were retrieved from the IUCN Red List of Threatened Species (https://www.iucnredlist.org/) and map was downloaded from https://www.diva-gis.org/gdata.

Table S2 :

 S2 List of study sites in Laos with coordinates

		SARS-CoV-2	SARS-CoV-2 non-neutralizing
		neutralizing sera		sera
	Extended data 1.	sp1		45	45
	Province	District	sp2 Village	Locality	3 Latitude	Longitude
	Vientiane Province	Pipistrellus Site 1 : Fueng Scotomanes Tylonycteris	sp3 sp. Na Ang ornatus fulvida malayana Noon Hin Hae	BANAL-236 Wuhan Pha Tong 8 Tham Nam cave Tham Pha 14 3 Pha Nouk Kok	1 18° 32.879'N 101° 58.938'E 18° 32.914'N 101° 58.459'E 1 18° 34.018'N 101° 58.338'E 14 18° 31.743'N 101° 58.648'E
	Grand Total		sp. Pha Louang	5 335 Tham Pha Louang 165	2 132 18° 30.495'N 101° 57.421'E 13 645
		Site 2 : Meth	Had Yao	Pha Keo Phong Mon 18° 54.426'N 101° 58.789'E
					Pha Khee Chear (Keo Tao)	18° 54.898'N 101° 59.048'E
			Houay Pa Mak	Pha Mum	18° 53.717'N 101° 56.957'E
					Pha Mum, Keo Tao	18° 54.787'N 101° 58.292'E
					Pha Yao, Keo Tao	18° 55.386'N 101° 57.723'E
	Oudomxay Province	Site 3 : Xay	Chom Ong	Chom Ong Cave	20° 43.071'N 101° 45.811'E
					Houay Om Loung	20° 43.111'N 101° 46.577'E
					Houay Pak Veuk	20° 41.647'N 101° 47.179'E
					Pa Noo Khom	20° 44.143'N 101° 46.770'E
		Site 4 : Na Mor	Na Thong	Nam Ook Hoo Cave	20° 52.396'N 101° 46.973'E
						20° 51.907'N 101° 46.926'E
							1

Table S3 :

 S3 Results from pan-coronavirus RT nested PCR and Sanger sequencing

	Anal Swab ID	Bat Genus	Bat Species	Genus	Subgenus	Closest sequence	% blastN similarity	District	Village	Locality
	BANAL-20-27	Rhinolophus	pusillus	co-infection Alpha/Betacoronavirus	unclassified Sarbecovirus Decacovirus/	Rhinolophus bat SARS-CoV-2 coronavirus HKU32 /	80%	Feuang (Site 1) Nonhinhae Pha Nok Kok
	BANAL-20-51	Hipposideros	pomona	Alphacoronavirus	Decacovirus	Bat coronavirus HKU10	98%	Feuang (Site 1) Nonhinhae Pha Nok Kok
	BANAL-20-52	Rhinolophus	malayanus	Betacoronavirus	Sarbecovirus	SARS-CoV-2	99%	Feuang (Site 1) Nonhinhae Pha Nok Kok
	BANAL-20-103	Rhinolophus	pusillus	Betacoronavirus	Sarbecovirus	SARS-CoV-2	99%	Feuang (Site 1)	Na Ang	Tham Pha
	BANAL-20-115	Hipposideros	pomona	Alphacoronavirus	Decacovirus	Bat coronavirus HKU10	98%	Feuang (Site 1)	Pha Luang	Tham Pha Luang
	BANAL-20-116	Rhinolophus	malayanus	Betacoronavirus	Sarbecovirus	SARS-CoV-2	98%	Feuang (Site 1)	Pha Luang	Tham Pha Luang
	BANAL-20-178	Hipposideros	pomona	Alphacoronavirus	Decacovirus	Bat coronavirus HKU10	98%	Feuang (Site 1)	Pha Luang	Tham Pha Luang
	BANAL-20-191	Chaerephon	plicatus	Alphacoronavirus	Pedacovirus	Porcine epidemic diarrhea virus (PEDV)	93%	Feuang (Site 1)	Pha Luang	Tham Pha Luang
	BANAL-20-197	Chaerephon	plicatus	Alphacoronavirus	Pedacovirus	Porcine epidemic diarrhea virus (PEDV)	93%	Feuang (Site 1)	Pha Luang	Tham Pha Luang
	BANAL-20-212	Chaerephon	plicatus	Alphacoronavirus	Pedacovirus	Porcine epidemic diarrhea virus (PEDV)	93%	Feuang (Site 1)	Pha Luang	Tham Pha Luang
	BANAL-20-213	Chaerephon	plicatus	Alphacoronavirus	Pedacovirus	Porcine epidemic diarrhea virus (PEDV)	93%	Feuang (Site 1)	Pha Luang	Tham Pha Luang
	BANAL-20-236	Rhinolophus	marshalli	Betacoronavirus	Sarbecovirus	SARS-CoV-2	99%	Feuang (Site 1)	Na Ang	Pha Tong
	BANAL-20-242	Rhinolophus	malayanus	Betacoronavirus	Sarbecovirus	SARS-CoV-2	99%	Feuang (Site 1)	Na Ang	Pha Tong
	BANAL-20-247	Rhinolophus	malayanus	Betacoronavirus	Sarbecovirus	SARS-CoV-2	99%	Feuang (Site 1)	Na Ang	Pha Tong

Table S5 .

 S5 Characteristics of the recombination breakpoints

	Breakpoint no	Downstream fragment no	Coordinates of the fragment	Breakpoint located in coding region?	Breakpoint position in the codon	Sequences before (aa)	Breakpoint sequence	Sequences after (aa)
	1		1 -1,820	yes (ORF1a)	3rd	GGTGGTGTTGTTCA (GGVVQ)	G	TTGACTTCGCAGTGG (LTSQW)
	2		1,821 -2,689	yes (ORF1a)	2nd	CCACTGGGCATTG (PLGID)	A	TTAGATGAGTGGAGT (LDEWS)
	3		3,971 -7,176	yes (ORF1a)	2nd	GTGGTTATACCTA (VVIPT)	C	TAAAAAGGCTGGTGGC (KKAGG)
	4		7,177 -11,461	yes (ORF1a)	1st	ATGTGGTTAATAATT (MWLIII)	A	ATCTTGTACAAATG (NLVQM)
	5		11,462 -12562	yes (ORF1a)	2nd	TTTGGCCTCTTTT (FGLFC)	G	TTTACTCAACCGCTAC (LLNRY)
	6		12,563 -17694	yes (ORF1a)	2nd	ACTGCTTGCACTG (TACTD)	A	TGACAATGCGTTAGCT (DNALA)
	7		17,695 -19,791	yes (ORF1ab)	2nd	ACCACTGAAACAG (TTETA)	C	AGCTCACTCTTGTAAT (AHSCN)
	8		19,792 -21,204	yes (ORF1ab)	2nd	CTCACTGTCTTTT (LTVFF)	T	TGATGGTAGAGTTGAT (DGRVD)
	9	10	21,205 -22,388	yes (ORF1ab)	2nd	TTATTTGACATGA (LFDMS)	G	TAAATTTCCCCTTAAA (KFPLK)
	10	11	22,389 -24,229	yes (S)	3rd	ACAGAATCTATTGT (TESIV)	T	AGATTTCCTAATATT (RFPNI)
	11	12	24,230 -27,343	yes (S)	2nd	CAAGACTCACTTT (QDSLS)	C	TTCCACAGCAAGTGCA (STASA)
	12	13	27,344 -27,800	yes (ORF7a)	2nd	GAGTGTGTTAGAG (ECVRG)	G	TACAACAGTACTTTTA (TTVLL)
	13	14	27,801 -28,463	yes (ORF8)	1st	CTTGTTTTCTTAGGA (LVFLG)	A	ATCATCACAACTGTA (IITTV)
	14	15	28,464 -29,597	yes (N)	3rd	AGAGCTACCAGACG (RATRR)	A	ATTCGTGGTGGTGAC (IRGGD)

Table S6 .

 S6 List of the MD simulations performed. The 3 production simulations of the SARS-CoV-2 RBD/hACE2 complex were initiated from the same equilibrated model with different random velocities. Production simulations of the BANAL-236 and BANAL-52/103 RBD/hACE2 complexes were initiated from different equilibrated models obtained, for each construct, from the top 3 scoring homology models.

Table S7 . Data collection and refinement statistics of the crystal structure of RBD-236 in complex with hACE2.

 S7

		Crystal 1
	Data collection
	Space group	P 65 2 2
	Cell dimensions
	a, b, c (Å)	341.53 341.53 68.15
	 ()	90 90 120
	Resolution (Å)	48.63 -2.94
		(3.047 -2.94) *
	Rsym or Rmerge	0.05642 (0.4364)
	I / I	11.42 (1.50)
	Completeness (%)	86.98 (12.53)
	Redundancy	2.0 (2.0)
	Refinement
	Resolution (Å)	48.63 -2.94
		(3.047 -2.94) *
	No. reflections
	Rwork / Rfree	0.214/0.227
	No. atoms	6610
	Protein	6418
	Ligand/ion	184
	Water	8
	B-factors	64.68
	Protein	63.72
	Ligand/ion	99.31
	Water	42.82
	R.m.s. deviations
	Bond lengths (Å)	0.004
	Bond angles ()	0.65

*Values in parentheses are for highest-resolution shell.

Acknowledgments:

The authors want to thank Souand Mohamed Ali, Nicolas Da Rocha, Angela Brisebarre, Thonglakhone Xaybounsou, and Somsanith Chonephetsarath for their help at the bench. Next-generation sequencing was performed with the help of Biomics Platform, C2RT, Institut Pasteur, Paris, France, supported by France Génomique (ANR-10-INBS-09-09), IBISA, and the Illumina COVID-19 Projects' offer. The work was granted access to the HPC resources of IDRIS under the allocation 2020-101592 made by GENCI. We thank the Ministry of Health and the Ministry of Natural Resources and Environments, Lao PDR, for their authorization of the field work and the Faculty of Environmental Science for its authorization of the field research collaboration. The work was funded by an Institut Pasteur "Covid Taskforce" and in part by the H2020 project 101003589 (RECOVER) and Labex IBEID (ANR-10-LABX62-IBEID) grants. Field and laboratory work at IP-Laos was also funded by a UK embassy grant (Grant No. INT 2021/LOV C19 02) and Luxembourg Development special grant (Grant No. LAO/030202324).

Data availability statement:

Sequences data that support the findings of this study have been deposited in the GenBank and GISAID databases with the following accession numbers: MZ937000 / EPI_ISL_4302644 (BANAL-52), MZ937001 / EPI_ISL_4302645 (BANAL-103), MZ937002 / EPI_ISL_4302646 (BANAL-116), MZ937003 / EPI_ISL_4302647 (BANAL-236), and MZ937004 / EPI_ISL_4302648 (BANAL-247). Raw sarbecovirus NGS reads and Sanger sequencing .ab1 files were deposited onto the NCBI SRA repository under the BioProjet PRJNA796968 (SAMN24959173-77). The crystal structure presented in this manuscript was deposited in the Protein Data Bank (PDB) with accession code 7PKI. The GROMACS topology and input files as well as the analysis scripts used are freely available on PLUMED-NEST (www.plumed-nest.org) under accession ID plumID:21.037.

Code availability statement:

MODELLER v. 10.1 (https://salilab.org/modeller/) and GROMACS v. 2020. 4 (http://www.gromacs.org) were used for MD simulations. XDS (v. feb 5 2021), CCP4 (v 7.0), COOT (v. 0.8.6) and PHENIX (v. 1.19.2-4158) softwares were used for X-ray diffraction data processing, model building and refinement of BANAL 236 RBD-hACE2 complex . PyMOL (v 2.4.2) was used for structural image rendering. PLUMED v. 2.7 (https://www.plumed.org); ROSETTA v. 3.11 (https://www.rosettacommons.org/); FoldX v. 4 (http://foldxsuite.crg.eu/); MDAnalysis v. 1.0.0 (https://www.mdanalysis.org/); MDTraj v. 1.9.5 (https://www.mdtraj.org/1.9.5/index.html); and in-house scripts available at https://github.com/maxbonomi/bat-MD were used for MD simulations.

hyphy 2.5.31; PhyML & MAFFT implemented through NGPphylogeny (https://ngphylogeny.fr/); Simplot 3.5.1 were used for phylogenetic and recombination analyses. MSSPE 1 was used for Betacoronavirus primer enrichment design and is available at https://github.com/chiulab/MSSPE-design.

MICROSEEK is an in-house pipeline that uses ALIENTRIMMER v. 2.0 (https://gitlab.pasteur.fr/GIPhy/AlienTrimmer) for read trimming / clipping, BBNORM from BBMAP v. 38.86 package (https://sourceforge.net/projects/bbmap/) for coverage normalization, MEGAHIT v. 1.2.9 (https://github.com/voutcn/megahit) for assembly, a in-house ORF-finder (https://figshare.com/articles/code/translateReads_py/7588592), and then DIAMOND v. 2.0.4 (https://github.com/bbuchfink/diamond/) and NCBI BLAST v. 2.12.0+ (ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.12.0/) both for sequence searching.

Author contributions

ST, KV, PTB, ME conceived the study design.

KV, BD, KL, NPDS, VX, PPa carried out sample procurement and bat species identification. SM, FD performed virus isolation, entry, and neutralization assays. BR, DC, PP prepared the NGS libraries and carried out the next-generation sequencing.

ST, BR, TB carried out the genome assembly, recombination, and phylogenetic analyses.

VL, SS, KL, NP implemented the pan-coronavirus PCR testing.

EBS and FR performed and analyzed the structure and binding studies. MB, YK, MN performed and analyzed the molecular dynamic simulations. ME and ST wrote the manuscript with inputs from all other authors.

Competing interests

All authors of the manuscript are cited as inventors of a U.S. Provisional patent application number 63/240.535 filed on September 03, 2021 entitled: NEW ISOLATED BAT SARS-COV ABLE TO INFECT HUMAN CELLS AND HARBORING RECEPTOR-BINDING DOMAINS CLOSE TO SARS-COV-2. The patent covers medical applications (diagnostic, vaccination) of the sequences described in the manuscript, more specifically in the part " Diversity of bat and coronavirus species" 0 20 40

Supplementary Methods

Additional analysis of the MD simulations of the BANAL-52/103 RBD/hACE2 complex. One of our 3 simulations of the BANAL-52/103 RBD/hACE2 complex (BANAL-52/103-CoV.1 in Extended data 5) showed large fluctuations of the RBD and RBD-hACE2 interface RMSDs. As a consequence of this larger conformational variability, the cluster analysis of BANAL-52/103-CoV.1 identified 2 major populated clusters (Extended data 5E). To exclude the possibility of having initialized the simulation from a poorly equilibrated structure, we performed a 10 times longer equilibration phase of the initial model followed by an additional 1 µs-long production simulation (BANAL-52/103-CoV.1*). The time series of the RBD-hACE2 interface RMSD in this additional simulation (Extended data 8, pink line) still displayed large RMSD fluctuations, especially towards the end of the simulation. After careful inspection of the trajectory, we concluded that the cause of the observed RMSD behavior in both BANAL-52/103-CoV.1 and BANAL-52/103-CoV.1* was the structural variability of the loop formed by RBD residues S443 to Y449 Extended data 8, insets, yellow residues). According to our definition, these residues were considered at the RBD/hACE2 interface; however, only a transiently formed hydrogen bond maintains this loop close to hACE2 (Extended data 7A). The removal of this stretch of 7 residues (out of 110 interfacial residues) from the calculation of the interface RMSD (Extended data 8, black and gray lines) clearly shows that the motion of the S443-Y449 loop was the cause of the large RMSD fluctuations observed in the BANAL-52/103-CoV.1 simulation. Amino-acid (lower) and nucleotide (upper) identity matrices of Laotian and representative human, bat, and pangolin sarbecoviruses. Individual genes and proteins were aligned with MAFFT, and identity matrices were constructed using CLC Main Workbench 21.0.4 (Qiagen). Matrices were colored according to the identity scale, from 25% (red) to 100% (green) of nucleotide or amino-acid identity. 20

Supplementary Tables