Arnaud Fichant

Arnaud Felten

Armel Gallet

Olivier Firmesse

Mathilde Bonis
email: mathilde.bonis@anses.fr

Identification of genetic markers for the discrimination of Bacillus thuringiensis within the Bacillus

Keywords: Bacillus thuringiensis, Bacillus cereus, foodborne outbreak, genome-wide association study

à 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

Bacillus cereus sensu lato, or the Bacillus cereus group (Bc), is composed of Gram-positive, facultative anaerobic, spore-forming, ubiquitous bacteria. This group comprises at least nine species: B. anthracis, B. cereus sensu stricto (Bc s.s.), B. thuringiensis (Bt), B. mycoides, B. pseudomycoides, B. weihenstephanensis, B. cytotoxicus, B. toyonensis and B. wiedmannii [START_REF] Miller | Bacillus wiedmannii sp. nov., a psychrotolerant and cytotoxic Bacillus cereus group species isolated from dairy foods and dairy environments[END_REF]. Furthermore, Bc members have very similar characteristics and highly conserved genomes (Ehling-Schulz et al., 2019). As a result, the taxonomic classification within the group is regularly revised, particularly due to the widespread use of whole-genome sequencing (WGS) which has led to higher resolution for the definition of new species [START_REF] Carroll | Keeping up with the Bacillus cereus group: taxonomy through the genomics era and beyond[END_REF][START_REF] Liu | Proposal of nine novel species of the Bacillus cereus group[END_REF].

Bc possesses an extensive number of virulence genes and is well known to be an agent responsible for foodborne outbreaks (FBOs). According to the European Food Safety Authority (EFSA), in 2019, Bc represented the leading cause of FBOs due to bacterial toxins in Europe [START_REF] Efsa | The European Union One Health 2019 Zoonoses Report[END_REF]. Bc can cause gastrointestinal disorders, such as diarrhea, vomiting, abdominal pain, or a combination thereof. The diarrhea is due to the ingestion of bacteria that produce enterotoxins, including three major pore-forming toxins: the non-hemolytic enterotoxin (Nhe) is found in approximately 85% to 100% of Bc strains, and hemolysin BL (Hbl) and cytotoxin K (CytK) are found in approximately 40 to 70% [START_REF] Dietrich | The Food Poisoning Toxins of Bacillus cereus[END_REF]. In contrast, the cereulide toxin, produced in food by some vegetative Bc, can lead to emetic syndromes, which predominantly cause vomiting. The ces operon coding for the emetic toxin synthesis genes is located on a megaplasmid in a restricted range of Bc s.s. strains [START_REF] Ceuppens | Regulation of toxin production by Bacillus cereus and its food safety implications[END_REF]. More rarely, Bc can also be responsible for extradigestive pathologies such as ocular and nosocomial infections [START_REF] David | Bacillus cereus endophthalmitis[END_REF][START_REF] Glasset | Bacillus cereus, a serious cause of nosocomial infections: Epidemiologic and genetic survey[END_REF].

Unlike other Bc bacteria, during sporulation, Bt strains have the ability to produce parasporal crystals containing toxins, some of which are extremely toxic toward insect larvae; Bt strains are therefore widely used as biopesticides in organic and conventional agriculture [START_REF] Roh | Bacillus thuringiensis as a specific, safe, and effective tool for insect pest control[END_REF]. These crystals are composed of Cry and Cyt protoxins (also called δ-endotoxins) for which several hundred different haplotypes have been described so far [START_REF] Crickmore | Revision of the nomenclature for the Bacillus thuringiensis pesticidal crystal proteins[END_REF]. Their classification has been recently revised based on protein sequence identity [START_REF] Crickmore | A structure-based nomenclature for Bacillus thuringiensis and other bacteria-derived pesticidal proteins[END_REF]. Upon ingestion of spores and crystals by larvae, larval digestive enzymes allow the release and activation of toxins from the crystals. Cry/Cyt toxins bind to specific receptors of the host midgut, forming pores in the cytoplasmic membrane of enterocytes, thereby triggering their death and creating breaches in the intestinal lining. In parallel, the favorable environment of the midgut supports the germination of spores and their entry in the internal milieu [START_REF] Palma | Bacillus thuringiensis toxins: an overview of their biocidal activity[END_REF]. Vegetative bacteria can proliferate, ultimately leading to septicemia and the death of the host. Four subspecies of Bt are commonly used in commercial products. Bt ssp. aizawai (Bta) and Bt ssp. kurstaki (Btk) and target Lepidoptera, Bt ssp. morrisoni (Btm) targets Coleoptera and Bt ssp. israelensis (Bti) is used against mosquitoes. Due to the ability of Bt to produce a large spectrum of insecticidal molecules, it is now considered as the leading microbial pesticide worldwide [START_REF] Lacey | Insect pathogens as biological control agents: back to the future[END_REF].

Interestingly, Bc and Bt share some virulence genes, particularly those encoding enterotoxins [START_REF] Dietrich | The Food Poisoning Toxins of Bacillus cereus[END_REF]. In addition, bacterial spores, in particular Bt spores, have the ability to persist in the environment and have been found on vegetables that have been treated with commercial products [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF][START_REF] Frederiksen | Occurrence of natural Bacillus thuringiensis contaminants and residues of Bacillus thuringiensis-based insecticides on fresh fruits and vegetables[END_REF][START_REF] Frentzel | Indications of biopesticidal Bacillus thuringiensis strains in bell pepper and tomato[END_REF]. Because Bc and Bt are genetically very closely related, there is no accurate test to distinguish between different Bc group members, meaning that discrimination between Bc s.s. and Bt is not straightforward. Nevertheless, optical microscopy can be used to detect parasporal insecticidal crystals (NF EN ISO 7932/Amd1), but this method requires expertise and several days to obtain analysis results. For instance, retrospective analyses in Canada [START_REF] Mcintyre | Identification of Bacillus cereus group species associated with food poisoning outbreaks in British Columbia, Canada[END_REF] and France [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF] identified Bt in FBOs that were initially attributed to Bc. In a report dealing with the possible involvement of Bt in foodborne infection events, EFSA called for the development of a simple method to differentiate Bt from other Bc members [START_REF] Efsa | Risks for public health related to the presence of Bacillus cereus and other Bacillus spp. including Bacillus thuringiensis in foodstuffs[END_REF], and a fortiori the Bt strains used as pesticides. In addition, the issue of food monitoring and, in particular, the acceptable dose of Bc/Bt in food is regularly examined by health authorities. The development of new tools for the specific identification of Bt strains would facilitate and extend this monitoring, with the aim of preventing and limiting the emergence of new intoxication cases.

Genetic proximity between Bt and other members of the Bc group, especially Bc s.s. [START_REF] Carroll | Keeping up with the Bacillus cereus group: taxonomy through the genomics era and beyond[END_REF]Ehling-Schulz et al., 2019), makes it difficult to develop identification methods; however, some genomic approaches may offer a solution. With the development of highthroughput sequencing (HTS), the rapid acquisition of complete genomes has been greatly facilitated. Based on pangenome analyses, recent studies have split the Bc group into three to five major clades [START_REF] Bazinet | Pan-genome and phylogeny of Bacillus cereus sensu lato[END_REF]Ehling-Schulz et al., 2019), whereas a previous classification defined seven phylogenetic groups, named I to VII and associated with different levels of psychrotolerance [START_REF] Guinebretière | Ecological diversification in the Bacillus cereus group[END_REF]. According to the literature, the entire Bt species is spread across several of these clades or groups, but most Bt strains of interest with regard to food safety (e.g. isolated from pesticides or FBOs) show less divergence [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF] and cluster with other strains in phylogenetic group IV, clade 2. In an attempt to distinguish them, four proximity clusters (named a to d) have been defined using an SNP calling approach.

The genome-wide association study (GWAS) approach, originally developed for human studies, has been adapted for the analysis of microorganisms [START_REF] Power | Microbial genome-wide association studies: lessons from human GWAS[END_REF]. The GWAS approach can establish a statistical link between a feature polymorphism (named hereafter "trait") and a molecular marker, such as a gene, a k-mer or a single nucleotide polymorphism (SNP). Several studies have already been carried out to correlate the presence of associated markers with a given phenotype, for example in the context of hostpathogen interactions [START_REF] Vila Nova | Genetic and metabolic signatures of Salmonella enterica subsp. enterica associated with animal sources at the pangenomic scale[END_REF], persistence and resistance to antibiotics [START_REF] Desjardins | Genomic and functional analyses of Mycobacterium tuberculosis strains implicate ald in D-cycloserine resistance[END_REF] as well as low-temperature growth phenotypes in food production [START_REF] Fritsch | Insights from genome-wide approaches to identify variants associated to phenotypes at pan-genome scale: Application to L. monocytogenes' ability to grow in cold conditions[END_REF]. Nevertheless, to date, no studies have identified markers specific to bacterial species or subspecies for the investigation on FBO etiological agents.

In this study, we searched for specific marker sequences for the identification of Bt and some Bt groups of interest, among other Bc. We performed a pan-GWAS analysis on 230 complete genomes. A total of 249 markers strongly associated with various traits (Bt species, Bt subspecies and Bt clusters a to d) were identified. Among them, 128 showed potential informativeness as genomic markers after a TBLASTN validation step. For each designated trait, one marker or a combination of up to six markers could identify all the Bt in the studied dataset. We wrote a Python script using 14 markers to allow accurate and automated identification of Bt from a genomic assembly.

Materials and methods

Whole-genome sequencing

The genomic DNA of 78 Bt/Bc isolates was extracted and sequenced as described in [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF]. Briefly, the KingFisher Cell and Tissue DNA kit (ThermoFisher) was used to isolate genomic DNA. DNA purity and concentration were determined using a Nanodrop Spectrophometer and a Qubit fluorimeter, respectively. Global DNA integrity was visualized on a 0.8% agarose gel (Seakem GTG™ Agarose) after migration for 2 h at 90 V. Library preparation was carried out using the Nextera XT DNA Library Prep kit (Illumina) and 150 bp paired-end sequencing of isolated DNA was performed by the Institut du Cerveau et de la Moëlle epinière, using a Nextseq500 sequencing system (Illumina). An in-house workflow called ARtWORK v1.0 [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF][START_REF] Mahamat Abdelrahim | Large-Scale Genomic Analyses and Toxinotyping of Clostridium perfringens Implicated in Foodborne Outbreaks in France[END_REF] was used to assemble reads with default parameters. The paired-end reads of the isolates are available in the PRJNA781790 BioProject on NCBI and the associated accession numbers are listed in Supplementary Table S1.

Dataset definition

A total of 286 genomes were used, corresponding to 215 FBO-associated Bc isolates, 18 commercial Bt strains [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF], and 53 Bc isolates of various origins, whose complete genomes were retrieved from the NCBI RefSeq database. In addition to the detection of parasporal insecticidal crystals under a microscope (NF EN ISO 7932/Amd1), Bt membership was confirmed by the presence of insecticidal protein-coding genes using the BtToxins_Digger pipeline v1.0.5 [START_REF] Liu | BtToxin_Digger: a comprehensive and high-throughput pipeline for mining toxin protein genes from Bacillus thuringiensis[END_REF]. Using R v4.0.1, a random draw was carried out to split the whole dataset into a study dataset (SD) and a validation dataset (VD) (Table S1), representing respectively 80% (n = 230) and 20% (n = 56) of the whole dataset.

The SD was used as an input for the GWAS analysis, and the VD was used to validate the selected genetic markers. The SD included the genomes of 144 strains belonging to Bt species, especially from the subspecies aizawai (n = 56) and kurstaki (n = 57), and 86 genomes belonging to other Bc, including at least one assembly for each of the nine distinct representative species of the group (Table 1) and for each of the seven previously defined phylogenetic groups [START_REF] Guinebretière | Ecological diversification in the Bacillus cereus group[END_REF] given in Table S1.

Pangenome analysis

Pangenome construction

An analysis was carried out to deduce the pangenome (i.e. core and accessory genes) from all 230 SD genomes (Fig. 1). The whole genome annotation of the assembly was performed using Prokka v1.13.3 [START_REF] Seemann | Prokka: rapid prokaryotic genome annotation[END_REF] with default parameters and GFF3 files were provided to Panaroo v1.2.3 (Tonkin-Hill et al., 2020) as input. Panaroo was run using the strict mode with default identity parameters and length difference thresholds (98%) for initial clustering of protein sequences with CD-HIT v4.8.1 [START_REF] Li | Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences[END_REF]. Then, close clusters were collapsed into putative families when they had at least 70% sequence identity.

The assignment of a gene cluster to a core genome was defined by the presence of the gene in at least 99% of the genomes (n = 228). As output, Panaroo produced a presence/absence matrix, used for the pan-GWAS analysis.

Genotype association

A three-level pan-GWAS analysis was performed to associate the presence or absence of a gene with a particular phenotype or trait. The first GWAS analysis (named L1 for Level 1) was conducted to search for genes specific to Bt species and the second analysis (named L2 for Level 2) for the Bt subspecies: Bt ssp. aizawai (Bta), Bt ssp. kurstaki (Btk), Bt ssp. israelensis (Bti) and Bt ssp. morrisoni (Btm). Then, specific genes for the SNP proximity clusters a, b, c and d within Bt subspecies (a, b for Bta and c, d for Btk) were investigated (analysis L3 for Level 3). These clusters were identified in a previous SNP study (see [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF]) and involve Bta or Btk strains that show strong genetic proximity with commercial strains used as insecticides. After targeting traits for each GWAS analysis, a gene search was carried out within the SD. Scoary v1.6.16 [START_REF] Brynildsrud | Rapid scoring of genes in microbial pan-genome-wide association studies with Scoary[END_REF] was run to identify statistically robust genes (Fisher's test) associated with each of the designated traits. For each gene, 1000 replicate permutations were conducted and p-values were adjusted by applying the Bonferroni correction [START_REF] Abdi | The Bonferonni and Šidák Corrections for Multiple Comparisons[END_REF] and the Benjamini-Hochberg procedure [START_REF] Wright | Adjusted P-Values for Simultaneous Inference[END_REF]. Only genes with p-values < 0.05, 100% specificity and at least 80% sensitivity for the screened traits were selected for the validation steps. In case of cluster-specific gene investigation, GWAS analysis (L3) was performed only on SD genomes with clear cluster attribution (n = 104), due to the impossibility to determine trait-specific genes when performed against all SD genomes. For each selected gene cluster, a protein sequence was extracted and tested as a candidate marker.

Core genome phylogeny

Core-genome sequences from Panaroo were aligned using MAFFT v7.471 [START_REF] Katoh | MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform[END_REF]. A phylogenetic tree of the SD was built with the maximum-likelihood method using iQtree v1.6.9 [START_REF] Nguyen | IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[END_REF] and the best substitution model GTR+F+R10 was determined with ModelFinder [START_REF] Kalyaanamoorthy | ModelFinder: fast model selection for accurate phylogenetic estimates[END_REF] on 286 DNA models. The core genome phylogeny was visualized using Phandango v1.3.0 [START_REF] Hadfield | Phandango: an interactive viewer for bacterial population genomics[END_REF], alongside pangenome metadata. Selected marker presence/absence information within the SD was added to the phylogenetic tree using iTol v6.1.2 [START_REF] Letunic | Interactive Tree Of Life (iTOL) v4: recent updates and new developments[END_REF].

In silico validation of genetic markers

To certify the reliability of the selected markers, we carried out a TBLASTN v2.7.1+ against the VD to validate candidate marker sequences meeting the required conditions. For each cluster associated with a given marker, a protein sequence was retrieved from the Panaroo output and screened for in the VD genomes. Markers were selected only when sequence coverage and identity were strictly higher than 90% compared with genomes displaying the same trait (intragroup), and when the sequence coverage and identity were strictly lower than 80% compared with genomes not associated with the specific trait (intergroup). This extra step made it possible to test the pan-GWAS selected markers in silico against a new dataset of genomes, and to confirm the sensitivity and the specificity of the markers. To enhance the identification of a genotype, a combination of validated makers associated with the best sensitivity was defined. Predicted genomic localization (i.e. chromosome vs. plasmid) of the gene encoding for the protein selected as a marker was carried out with 8 TBLASTN using protein sequences and a Bt plasmid database (n = 342) collected from PLSBD [START_REF] Galata | PLSDB: a resource of complete bacterial plasmids[END_REF].

Results

Bacillus cereus group pangenome

To identify specific Bt markers associated with different traits of interest using a pan-GWAS approach with Scoary, the pangenome of the Bc group was inferred using Panaroo on SD, resulting in a total of 39,021 clusters of genes. The core genome of the 230 Bc complete genomes was composed of 1854 genes (representing about 5% of all genes), leaving 37,167 accessory genes (Fig. 2). Among the accessory genes, 1695 were considered soft core genes (presence in 95 to 99% of all genomes), 4358 shell genes (presence in 15 to 95% of all genomes) and 31,114 cloud genes (presence in less than 15% of all the genomes). A rarefaction curve of the 230 genomes (Supplementary Fig. S1) suggests the existence of an open pangenome, as previously assumed for Bc [START_REF] Bazinet | Pan-genome and phylogeny of Bacillus cereus sensu lato[END_REF][START_REF] Zwick | Genomic characterization of the Bacillus cereus sensu lato species: backdrop to the evolution of Bacillus anthracis[END_REF]. The presence-absence matrix showed that the B. cytotoxicus population, which roots the tree, seemed to possess a reduced number of accessory genes. In contrast, the Bta and Btk subspecies, whose genomes represent highly conserved populations, displayed a greater number of accessory genes.

Pan-GWAS analysis

Gene-based GWAS

We used Scoary to associate specific genes with the presence or absence of specific traits with a GWAS analysis. Three GWAS analyses were run to identify a large number of candidate genes for each trait (Table 2). A total of 1304 genes were significantly associated with the traits of interest. Interestingly, the majority of the identified genes were associated with subspecies israelensis and morrisoni in the GWAS L2 analysis, with 809 and 246 genes, respectively. However, this large number of associated genes can be explained by the low number of genomes available for these two subspecies in the dataset reducing the analytical strength. Given the low diversity in the dataset for these subspecies, we did not pursue the search for specific markers for these subspecies.

Among the other 249 candidate genes, 131 genes met the specificity and sensitivity requirements for the L1 and L2 GWAS analyses, corresponding respectively to the Bt species trait and the Bt subspecies traits of the two remaining subspecies (i.e. Btk and Bta).

Accordingly, 99 genes were specific to the Bt species and 27 to the Bta subspecies.

Interestingly, only 5 genes were specific to Btk suggesting relatively high diversity within the Btk subspecies. For the L3 GWAS analysis performed on a smaller panel of genomes, 118 genes were identified among the four traits corresponding to clusters a, b, c and d with 10, 14, 46 and 48 genes, respectively.

Validation of GWAS results

We then validated the markers selected using GWAS by comparing them in silico to the VD using TBLASTN. For each selected candidate gene associated with a specific trait in the GWAS analysis, a protein sequence was extracted and compared to all VD genomes that shared or did not share the trait in question. Out of the 249 highly associated genes obtained with Scoary, 128 protein sequences met the criteria of coverage, specificity and sensitivity of the TBLASTN analysis (Table 2). All validated marker sequences had a sensitivity range between 80% and 100% on both SD and VD datasets. Then, when a single marker was not sufficient, combinations of markers were tested to identify all the associated genomes for each trait. Table 3 lists 14 markers associated with low p-values and whose combination could detect all genomes of a given trait (i.e. 100% sensitivity).

Functions and distribution of specific markers

To predict whether the genes coding for the markers were localized or not on a plasmid, a TBLASTN alignment on a Bt plasmid database was carried out (Table 3). When the function was known, the closest protein identified in Uniprot (The UniProt, 2017) was indicated. The Bt species can thus be identified with a combination of only six markers. For example, the cwlA gene located in an operon along with cwlB, encodes a peptidoglycan hydrolase, known for its role in cell lysis during Bt spore release [START_REF] Yang | Transcriptional regulation and characteristics of a novel N-acetylmuramoyl-L-alanine amidase gene involved in Bacillus thuringiensis mother cell lysis[END_REF]. Regarding the markers that emerged from the GWAS L2 analysis for Bt subspecies, all Bta genomes (n = 56) could be identified based on the apr gene, which encodes a specific chromosomal alkaline serine protease, also called subtilisin in Bacillus subtilis [START_REF] Park | Bacillus subtilis subtilisin gene (aprE) is expressed from a sigma A (sigma 43) promoter in vitro and in vivo[END_REF]. For Btk, two genes encoding hypothetical proteins (group_27293 and group_27336) were needed to identify the 57 genomes present in the SD. Moreover, both genes were located on plasmids suggesting the presence of a common plasmid in these subspecies. Regarding the cluster-specific markers identified with the GWAS L3 analysis, clusters a, c and d could each be identified with a single chromosomal marker, respectively group_10114 encoding a hypothetical protein, the repressor lexA, involved in the SOS system [START_REF] Au | Genetic composition of the Bacillus subtilis SOS system[END_REF] and the gene encoding a subunit of Clp protease [START_REF] Krüger | The clp proteases of Bacillus subtilis are directly involved in degradation of misfolded proteins[END_REF]. Cluster b, like Btk, needed two markers for the identification of all genomes: the rapF gene, involved in bacterial competence [START_REF] Bongiorni | Synergistic regulation of competence development in Bacillus subtilis by two Rap-Phr systems[END_REF], combined with a plasmid gene of unknown function (group_20667).

For all hypothetical protein-coding genes, we attempted to annotate them with the 3D structure using SWISS-MODEL [START_REF] Waterhouse | SWISS-MODEL: homology modelling of protein structures and complexes[END_REF]. However, it was not possible to assign new putative functions to them (Supplementary Table S3).

10 Phylogenetic reconstruction of the 1854 core genes from the 230 genomes (Fig. 3) showed a clear distinction between most Bc genomes, including Bt genomes, and five genomes of the B. cytotoxicus species, which is known to be the most distant group member (Ehling-Schulz et al., 2019). In contrast, the genomes of the two subspecies Btk and Bta were located on the outermost branch of the tree compared with B. cytotoxicus genomes. Unsurprisingly, other Bt subspecies or unknown subspecies genomes were scattered around the tree. Some of them, such as MTC28 (NCBI accession: NC_018693.1) and Bt finitimus YBT020 (NCBI accession: NC_017200.1) were phylogenetically very distant.

The 14 Bt-specific markers identified with the three levels of GWAS analysis and validated by TBLASTN were mapped on the phylogenetic tree (Fig. 3). Despite the great distance between some Bt strains, a combination of six markers can identify all of them. Similarly, with two markers specific to Btk and cluster b, all associated genomes can be retrieved. Only one marker was needed for the identification of Bta, clusters a, c and d.

Workflow for Bt identification

Based on the 14 selected markers, a workflow was developed to distinguish Bt from the genome of a putative Bc isolate (Fig. 4). We wrote an in-house script in Python 3 to automate the screening of the markers in genome assemblies. Using BLASTX, the script detects the markers for Bt, Bt subspecies and Bt clusters using the same thresholds as the TBLASTN validation (90% identity and coverage) in the user-input genomes. The output indicates the presence (+) or absence (-) of all trait-associated markers for each genome query. Then, according to the workflow pathway, based on the presence/absence of the markers within the genome, a putative Bt identification is given. The script is freely available at https://github.com/afelten-Anses/Bt_typing.

Discussion

To differentiate Bt from other members of the Bc group and especially in the context of food poisoning, we analyzed a large database of Bc to identify genetic markers specific to Bt. The construction of the Bc pangenome from 230 Bc genomes showed that it was composed of 39,021 protein-coding genes and that the core genome corresponded to 1854 genes. These results differ somewhat from those previously published for the Bc pangenome. A previous pangenome study of 114 Bc highlighted 59,989 genes, but only 600 were described as core genes [START_REF] Bazinet | Pan-genome and phylogeny of Bacillus cereus sensu lato[END_REF]. Normally, increasing the number of genomes in a dataset can lead to a decrease in the core genome size along with an increase in the accessory genome size.

However, the difference of tools used for pangenome deduction here can explain this large difference in gene numbers, as already highlighted in a genome-wide analysis of Clostridium difficile pangenome [START_REF] Knight | Major genetic discontinuity and novel toxigenic species in Clostridioides difficile taxonomy[END_REF]. For instance, the previous popular pangenome analysis software, named Roary [START_REF] Page | Roary: rapid large-scale prokaryote pan genome analysis[END_REF], uses a default threshold of 95% identity (versus 98% with Panaroo) to create clustered protein sequences. Nonetheless, Panaroo, unlike Roary, clusters potential close family genes together. Moreover, Roary does not take assembly and annotation errors into account, thereby resulting in an overestimation of the total gene number and thus a reduction of the core genome size, especially for datasets with high diversity as in the Bc group.

In addition, many Bt genomes, especially from the Btk and Bta subspecies, used for this study comprised a large number of accessory genes, which can be explained by the high diversity of the genomes in these subspecies in the dataset and/or the presence of mobile genetic elements, particularly plasmids. Several studies have already highlighted that the number of plasmids is high in Bt in comparison to Bc [START_REF] Fayad | Diversity of Bacillus cereus sensu lato mobilome[END_REF]. A comparative genomic analysis study has shown that Bt strains with high insecticidal potency harbor genes promoting infection, immune evasion and nutrient access that may play a key role in entomopathogenicity and host adaptation [START_REF] Liu | Genomic insights into the taxonomic status of the Bacillus cereus group[END_REF]. Toxin-carrying plasmids, which represent an important part of the Bt plasmid pool, have also been shown to be involved in cellular functions such as germination, sporulation and horizontal gene transfer [START_REF] Fayad | Diversity of Bacillus cereus sensu lato mobilome[END_REF][START_REF] Gillis | Role of plasmid plasticity and mobile genetic elements in the entomopathogen Bacillus thuringiensis serovar israelensis[END_REF]. For the Bti and Btm subspecies, for which only a limited number of genomes were included in this study, no gene validation step could be performed to define specific markers. Nevertheless, our analysis led to the identification of genes only present in a restricted population of Bti (18SBCL211A, 18SBCL484A, Bt_israelensis_4Q1, Bt_israelensis_HD789, Bt_israelensis_AM6552) or in one strain of Btm (Bt_morrisoni_BGSC 4AA1), whose specificity deserves further investigation.

A total of 249 candidate genes were selected based on their sensibility and specificity for a given trait. As recommended for GWAS analyses [START_REF] San | Current Affairs of Microbial Genome-Wide Association Studies: Approaches, Bottlenecks and Analytical Pitfalls[END_REF], the results obtained were validated by testing them on a test dataset (here, the VD). The fact that almost 50% of the selected markers did not pass this step demonstrates the importance of testing results.

To extend identification to other traits or Bt genomes, marker selection preferentially targeted genes with a chromosomal location and a known or predicted function. Of the 128 validated genes, only 21 had predicted functions, annotated with Prokka [START_REF] Seemann | Prokka: rapid prokaryotic genome annotation[END_REF] (Supplementary Table S3). At least one gene associated with a known function and located on the chromosome was identified for L1-and L2-associated traits (Bt species and subspecies, respectively), with the exception of Btk. However, we showed that Btk could be identified based on the combination of two plasmid markers, which underlines the importance of plasmid content, especially for the discrimination of closely related Bt strains.

12

The phylogenic visualization of the SD dataset illustrates the challenge -even when using genomic approaches -of differentiating Bt populations from each other or from other Bc group members, due to their close genetic proximity [START_REF] Bazinet | Pan-genome and phylogeny of Bacillus cereus sensu lato[END_REF][START_REF] Carroll | Keeping up with the Bacillus cereus group: taxonomy through the genomics era and beyond[END_REF]. For example, L3 analysis performed for cluster identification could not be performed on the whole SD. The impossibility of identifying specific genes revealed a limitation of gene-based GWAS when comparing extremely close groups. Fortunately, the sequential use of a combination of markers allows the use of cluster-specific markers after confirmation of Bt membership and then one of two subspecies, Bta or Btk. Furthermore, the use of a highly diverse dataset in this study allowed the identification of specific and high sensitivity markers. For example, among the Bt markers identified, the cwlA gene was found in 93.8% of SD Bt genomes (n = 144), making it a very reliable marker compared with the previous proposed detection system [START_REF] Chang | PCR assay of the groEL gene for detection and differentiation of Bacillus cereus group cells[END_REF][START_REF] Chelliah | Whole genome sequence of Bacillus thuringiensis ATCC 10792 and improved discrimination of Bacillus thuringiensis from Bacillus cereus group based on novel biomarkers[END_REF][START_REF] Dzieciol | A novel diagnostic real-time PCR assay for quantification and differentiation of emetic and non-emetic Bacillus cereus[END_REF]. Furthermore, this previous system did not include commercial Bt strains, unlike the analyses performed for this study.

Nevertheless, it is noteworthy that the identified markers refer to a specific dataset, meaning they may not allow the identification of the entire Bt species. In addition, we cannot exclude the existence of some non-sequenced strains, in particular divergent strains not included in our dataset that may not possess one or more of the selected markers. However, the workflow developed here can easily be adapted with the addition of new markers or modification of existing markers, to detect additional Bt strains of interest.

The 2016 EFSA report [START_REF] Efsa | Risks for public health related to the presence of Bacillus cereus and other Bacillus spp. including Bacillus thuringiensis in foodstuffs[END_REF] highlighted the need for the development of new methods to differentiate Bt from Bc s.s. Here, to differentiate Bt from other members of the Bc group, particularly in the context of FBOs, our work shows that a combination of six markers can identify Bt species, three for two major subspecies of interest in food safety, and five for the proximity clusters of pesticide strains frequently used in agriculture. Based on these 14 markers, we developed a new tool, associated with a workflow (Fig. 4) and a script, to predict the identity putative Bt among Bc isolates. Moreover, the gene-based GWAS approach developed here demonstrated that the four proximity clusters associated with commercial Bta or Btk strains likely belong to Bt [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF]. Currently, the workflow cannot be used to distinguish certain Bt strains (particularly the insecticide strains belonging to the kurstaki subspecies) within the same cluster. Gene-based approaches are limited in their specific identification of commercial Bt strains used in agriculture. New approaches, notably based on SNP calling, can be conducted to search for differentiation methods at the strain level. In this perspective, our workflow can be easily updated and optimized, if necessary.

With the routine use of sequencing methods in laboratories and the significant development of high-throughput sequencing techniques, the use of computational tools as a complementary method for the identification of bacterial species may prove to be a valuable asset, especially in the context of food poisoning. This complementary identification method can be used to quickly assign a Bc strain to the Bt species, as well as to challenge false positives and false negatives resulting from microscope searches for protein crystals, possibly due to ambiguous phenotypes, misinterpretation of results, or loss of plasmids carrying the cry genes. With different levels of analysis, possible assumptions on the origin of the isolates can be made. qPCR detection systems hold promise for the specific detection and quantification of Bt in food matrices, an important, as-yet-unattained goal for monitoring Bt in food.

Figure captions

Fig. 1. Flowchart of the pan-GWAS analysis. Pangenome construction was first performed using Panaroo [START_REF] Tonkin-Hill | Producing polished prokaryotic pangenomes with the Panaroo pipeline[END_REF] to obtain an estimation of the core and accessory genomes of the study dataset (SD). Then, the presence/absence matrix output from Panaroo was used for the GWAS analysis performed using Scoary, thereby providing a list of associated genes for each trait. Only genes with at least 80% sensitivity and 100% specificity for the SD were selected. A TBLASTN validation step of the selected genes was performed after recovering the protein sequences from Panaroo. The genes were validated when the respective sequence criteria of coverage and identity were met and when exhibiting at least 80% sensitivity and 100% specificity for the validation dataset (VD). The core-gene alignment obtained from Panaroo was used to construct the phylogeny of the 230 genomes, using IQtree [START_REF] Nguyen | IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[END_REF]. (SD) genomes. Visualization and annotation were performed using iTol [START_REF] Letunic | Interactive Tree Of Life (iTOL) v4: recent updates and new developments[END_REF]. Filled circles indicate the presence of the 14 Bt markers for the successive GWAS analysis levels (L1, L2 and L3) in the corresponding SD genomes. 3). A marker is considered present when its sequence coverage and its identity are both higher than 90%

in the assembly tested. The workflow was automated using a Python script available online (https://github.com/afelten-Anses/Bt_typing).

Fig. 2 .

 2 Fig. 2. Phylogenetic tree of 230 Bacillus cereus (Bc) genomes alongside the gene presence/absence variation matrix. The maximum-likelihood phylogeny (GTR+F+R10 substitution model) was performed with IQtree (Nguyen et al., 2015) using the Panaroo (Tonkin-Hill et al., 2020) core-gene alignment of the study dataset (SD) genomes and visualized with Phandango[START_REF] Hadfield | Phandango: an interactive viewer for bacterial population genomics[END_REF]. The presence-absence matrix obtained with Panaroo identified 1854 core genes (present in 99% of SD genomes) and 37,167 accessory genes for the 230 Bc genomes.

Fig. 3 .

 3 Fig. 3. Core genome phylogeny of 230 Bacillus cereus (Bc) genomes. The maximumlikelihood phylogeny (GTR+F+R10 substitution model) was performed with IQtree (Nguyen et al., 2015) using a Panaroo (Tonkin-Hill et al., 2020) core-gene alignment of the study dataset

Fig. 4 .

 4 Fig. 4. Workflow for Bacillus thuringiensis (Bt) identification. The proposed decision tree is used to determine if a putative Bacillus cereus (Bc) isolate belongs to the Bt species (L1), to the Bt ssp. kurstaki (Btk) or aizawai (Bta) subspecies (L2) and to the clusters a, b, c, and d (L3), based on the 14 marker protein sequences previously identified (see Table3). A marker

Table 3 .

 3 List of 14 genes with strong evidence and best combinations for Bacillus thuringiensis (Bt) trait identification. 611

	GWAS analysis	Gene name	Replicon	Annotation	UniprotKB	Trait	Marker Sensitivity (%)	Cumulative Sensitivity (%)	p-value	NCBI accession number
		cwlA	Chromosome	N-acetylmuramoyl-L-alanine amidase CwlA	P24808	Bt	93.8	93.8	1.72E-51	WP_021728236.1
	L1 L2	intQ group_3916 group_20749 Chromosome Chromosome Chromosome group_20361 Plasmid sdpR Plasmid apr Chromosome group_27293 Plasmid group_27336 Plasmid group_10114 Chromosome	Putative defective protein IntQ Hypothetical protein Hypothetical protein Hypothetical protein Transcriptional repressor SdpR Subtilisin Hypothetical protein Hypothetical protein Hypothetical protein	P76168 N/A N/A N/A O32242 P04189 N/A N/A N/A	Bt Bt Bt Bt Bt Bta Btk Btk Cluster a	91.7 88.9 86.8 84.7 84.0 100.0 98.2 89.5 100.0	94.4 95.1 97.2 98.6 100.0 N/A 98.2 100.0 N/A	1.64E-48 3.29E-47 6.36E-45 8.43E-43 4.00E-42 9.10E-53 5.83E-52 1.05E-43 4.48E-15	WP_000237488.1 WP_000858032.1 WP_042596929.1 WP_002101540.1 WP_000998670.1 WP_021728520.1 WP_003273526.1 WP_001293418.1 WP_000415284.1
		rapF	Chromosome	Response regulator aspartate phosphatase F	P71002	Cluster b	92.5	92.5	4.82E-25	WP_050062578.1
	L3	group_20667 Plasmid lexA Chromosome clpP1 Chromosome	Hypothetical protein LexA repressor ATP-dependent Clp subunit 1	N/A P31080 B0B803	Cluster b Cluster c Cluster d	82.5 100.0 100.0	100.0 N/A N/A	1.34E-20 4.31E-24 2.13E-26	WP_131256056.1 AHZ54004.1 WP_000791073.1

Acknowledgements

This project was funded and supported by ANSES (PPV project "OSABt") via the tax on sales of plant protection products. The proceeds of this tax are allocated to ANSES to finance the development of a system to monitor the adverse effects of plant protection products, called "phytopharmacovigilance" (PPV), established by the French Act on the Future of Agriculture, Food and Forestry of 13 October 2014. A. Fichant was funded by the Consortium Biocontrol (BACILLUS project).

We thank all the district veterinary and food analysis laboratories for carrying out Bc detection and transmitting isolates to the Food Safety Laboratory.

Tables

 Table 1.Composition of the study dataset (SD).

Supplementary data

Supplementary Fig. S1. Rarefaction curve of gene diversity as a function of the number of samples generated with non-parametric incidence-based estimator jack1 [START_REF] Burnham | Estimation of the size of a closed population when capture probabilities vary among animals[END_REF]. S1. List of the 286 genomes used in this study. The genome set was divided into two datasets, a study dataset (SD) (n=230) and a validation dataset (VD) (n=56).

Supplementary Table

a Clusters were determined by SNP calling, using iVARCall2 [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF]. b All the sequencing data used in this study are associated with two BioProjects: PRJNA547495 [START_REF] Bonis | Comparative phenotypic, genotypic and genomic analyses of Bacillus thuringiensis associated with foodborne outbreaks in France[END_REF] and PRJNA781790 (this study). The attribution to panC groups was performed after partial sequencing of panC gene, and according to Guinebretière et al., [START_REF] Guinebretière | Ecological diversification in the Bacillus cereus group[END_REF]. Abbreviations: FBO=foodborne outbreak. S2. Search for 3D structure prediction for workflow markers of Supplementary Table S3. List of the 21 annotated genes found with GWAS analysis and validated using TBLASTN.

Supplementary Table


Gwénaël Jan
email: gwenael.jan@inrae.fr

Florian Tarnaud

Fillipe Luiz

Rosa Do Carmo

Nassima Illikoud

Fanny Canon

Julien Jardin

Fanny Guyomarc'h

Valérie Briard-Bion

Valérie Gagnaire

Data from a proteomic comparative analysis highlight differential adaptation of Lactobacillus delbrueckii subsp. bulgaricus to cow milk versus to soy milk environments

Keywords: microbiology, dairy starter, yogurt starter, probiotic, stress adaptation

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

The article presents a proteomic dataset generated by a comparative analysis, using gel-free nanoLC-MS/MS, of the cellular proteome of Lactobacillus delbrueckii subsp. bulgaricus , a yogurt starter, when cultivated in soy milk versus in cow milk. The CIRM-BIA1592 strain was cultivated in the aqueous phase of soy milk, or of cow milk. Whole-cell proteins were extracted, trypsinolyzed and analyzed by nano LC-MS/MS, prior to identification and to classification by function using the X!Tandem pipeline software and the proteomic data from NCBI.nlm.nigh.gov. Quantification of the proteins was moreover performed to evidence changes in their expression, depending on the culture medium. Data are available via Pro-teomeXchange with the identifier PXD033905 ( http://www. proteomexchange.org/ ). This article is related to the research article entitled "The stressing life of Lactobacillus delbrueckii subsp. bulgaricus in soy milk", by G.Jan et al. in Food Microbiology, 2022. This proteomic differential analysis indeed revealed major modulation of the stress proteome, with many stress proteins upregulated in the soy environment.

Value of the Data

• These data constitute an important step in understanding the molecular mechanisms responsible for the adaptation of yogurt starters to the soy environment. • These data should be considered in the development of innovative fermented soy products. • These data evidence that the yogurt starter L. bulgaricus is stressed when cultivated in soymilk, compared to milk.

Data Description

This article presents a dataset generated during a comparison of the cellular proteome of the yogurt starter Lactobacillus delbrueckii subsp. bulgaricus cultivated in (1) cow milk environment and (2) soy milk environment. The list of all the cellular proteins that were identified in this study is provided in a table entitled "protein_quantification_results.csv" available via Pro-teomeXchange with the identifier PXD033905. Project Name: Data from a proteomic comparative analysis highlight differential adaptation of Lactobacillus delbrueckii subsp. bulgaricus to cow milk versus to soy milk environments. Project accession: PXD033905.

The abundance ratios of the differentially expressed proteins, displaying a significantly different abundance when bacteria were cultivated in soy versus in cow milk environment, classified in different COG categories, are published in the following article: [START_REF] Jan | The stressing life of Lactobacillus delbrueckii subsp. bulgaricus in soy milk[END_REF] . The stressing life of Lactobacillus delbrueckii subsp. bulgaricus in soy milk. Food Microbiol 106:104042. https: //doi.org/10.1016/j.fm.2022.104042 .

For each biological and each technical replicate, a raw data file is presented in the deposited dataset as mentioned above.

This article presents data from the research article entitled "The stressing life of Lactobacillus delbrueckii subsp. bulgaricus in soy milk". The cellular proteome, according to the growth medium, reveals major induction of stress proteins, in the soy environment.

Experimental Design, Materials and Methods

Strain and Pre-Cultures

De Man, Rogosa, Sharpe (MRS) medium [START_REF] De Man | A medium for the cultivation of lactobacilli[END_REF] was used to grow liquid precultures of Lactobacillus delbrueckii subsp. bulgaricus CIRM-BIA1592 at 42 °C without agitation. This strain was provided by CIRM-BIA, which is an International Centre for Microbial Resources dedicated to bacteria of food interest called CIRM-BIA, (INRAE Rennes, France, https://collection-cirmbia.fr/ ). These precultures were used to inoculate cultures, either in MRS, or in ultrafiltrate media, as described below.

Culture in Cow Milk or Soy Milk Ultrafiltrate

Raw cow milk was skimmed and then ultrafiltrated on the INRAE STLO dairy platform as described previously [ 3 , 4 ]. In that aim, we used a filtration pilot, containing a ceramic membrane. This membrane has a molecular weight cut-off point of 8 kDa. This ultrafiltration concentrates dairy proteins in the retentate and generates a permeate, MUF (milk ultrafiltrate), which constitutes the aqueous phase of cow milk. The overall composition of MUF was as follows: carbohydrate 5%; non-protein nitrogen 0.28%, minerals 0.75% and dry matter 6.14%. As an additional nitrogen source, 5 g.L -1 food grade casein hydrolysate (Casein Peptone Plus, Organotechnie, La Courneuve, France) was added. After adjustment of pH to 7.0 using NaOH, the MUF medium was sterilized using 0.2 μm filters (Nalgene, Roskilde, Denmark) prior to storage at 4 °C.

Soy milk ultrafiltrate (SUF) is a by-product of industrial concentration of soy proteins during the manufacture of soymilk. It was provided by a local company, Sojasun Technologies Triballat Noyal (Noyal-sur-Vilaine, France) and prepared according to patent N °EP 1 983 844 B1 [START_REF] Efstathiou | Method for production of active extracts from soya beans and corresponding uses of the produced extracts[END_REF] . This SUF is commercially available under the designation BASOSOY (Triballat ingredients Triballat, Noyal-sur-Vilaine, France). The overall composition of SUF was as follows: carbohydrates 2.5%, protein 0.55%, non-protein nitrogen 0.20%, mineral 1%, dry matter 5%. As an additional nitrogen source, 5 g.L -1 food grade soy hydrolysate (Bacto Soytone, BD Bioscience) was added. After adjustment of pH to 7.0 using NaOH, the SUF medium was streilized by autoclaving (110 °C, 10 min) prior to sterile centrifugation (28,0 0 0 × g, 30 min) and filtration in order to remove insoluble compounds. It was then sterilized using 0.2 μm filters (Nalgene, Roskilde, Denmark) prior to storage at 4 °C.

In this work, 4 independent cultures of Lactobacillus delbrueckii subsp. bulgaricus CIRM-BIA1592 were prepared at 42 °C without agitation in MUF. In addition, 4 independent cultures of Lactobacillus delbrueckii subsp. bulgaricus CIRM-BIA1592 were prepared at 42 °C without agitation in SUF. Cultures were incubated until the end of growth and of acidification, i.e. during 16 h. These 8 cultures were then subjected to whole-cell label-free proteomic analysis as described below. For each culture (4 biological replicates), 3 mass spectrometry analyses (3 technical replicates) were performed.

Label-Free Proteomics

Label-free proteomics was performed as described in details in the corresponding research article [START_REF] Jan | The stressing life of Lactobacillus delbrueckii subsp. bulgaricus in soy milk[END_REF] and adapted from a previous study [START_REF] Gaucher | Benefits and drawbacks of osmotic adjustment in Propionibacterium freudenreichii[END_REF] . Briefly, lactobacilli were grown either in milk or soy ultrafiltrate and then lysed in lysis solution (50 mM Tris-HCl [pH 7.5], 0.3% SDS, 200 mM dithiothreitol (DTT), 0.4 mM phenyl methyl sulfonyl fluoride, PMSF). Protein extracts were cleaned and quantified using the 2-D Clean-Up kit (GE Healthcare) and the 2-D Quant Kit, respectively. Whole-cell tryptic digestion was done using Sequencing Grade Modified Trypsin (Promega, Madison, USA) as described previously [START_REF] Huang | Growth in hyper-concentrated sweet whey triggers multi stress tolerance and spray drying survival in Lactobacillus casei BL23: from the molecular basis to new perspectives for sustainable probiotic production[END_REF] prior to addition of spectrophotometricgrade trifluoroacetic acid (TFA) (Sigma-Aldrich, USA).

Nano-LC-MS/MS was done as previously described [START_REF] Huang | Growth in hyper-concentrated sweet whey triggers multi stress tolerance and spray drying survival in Lactobacillus casei BL23: from the molecular basis to new perspectives for sustainable probiotic production[END_REF] using a nano RSLC Dionex U30 0 0 system fitted to a Q-Exactive mass spectrometer (Thermo Scientific, San Jose, USA). Spectra of eluted peptides were recorded in full MS mode, selected between 250-20 0 0 m/z, with a resolution of 70,0 0 0 at m/z 200.

Identification of proteins was as previously described [START_REF] Huang | Growth in hyper-concentrated sweet whey triggers multi stress tolerance and spray drying survival in Lactobacillus casei BL23: from the molecular basis to new perspectives for sustainable probiotic production[END_REF] , peptides being identified using the X!Tandem pipeline software (Langella et al., 2017), querying the proteome of L. delbrueckii subsp. bulgaricus (proteome UP0 0 0 0 01259 downloaded from uniprot.org).

Protein quantification was as described [START_REF] Huang | Growth in hyper-concentrated sweet whey triggers multi stress tolerance and spray drying survival in Lactobacillus casei BL23: from the molecular basis to new perspectives for sustainable probiotic production[END_REF] , each peptide being quantified using the Mass-ChroQ software prior to data treatment and statistical analysis under the R software (R 3.2.2, Project for statistical computing) using the specific R package called 'MassChroqR'. Variance analysis was performed on proteins with a minimum peak ratio of 1.5 for peak counting analysis. We considered significantly different proteins with an adjusted p -value < 0.05.

Proteins with a significant ( p < 0.05, ANOVA) change in expression of ≥ 2-fold (log2 ratio ≥ 1.5) were considered differentially expressed. We generated a volcano plot to illustrate differentially expressed proteins in L. delbrueckii subsp. bulgaricus . Functional annotation and Clusters of Orthologous Groups (COGs) were done with eggNOG-mapper v2 web tool (Huerta-Cepas et al., 2017, 2019).

The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE [START_REF] Perez-Riverol | The PRIDE database resources in 2022: a hub for mass spectrometry-based proteomics evidences[END_REF] partner repository with the dataset identifier PXD033905.

CRediT Author Statement

  Jan: Writing -review & editing, Visualization, Methodology, Supervision; Florian Tarnaud: Writing -review & editing, Methodology; Fillipe Luiz Rosa do Carmo: Writing -review & editing, Data curation; Nassima Illikoud: Writing -review & editing, Data curation; Fanny Canon: Writing -review & editing, Methodology; Julien Jardin: Writing -review & editing, Methodology; Valérie Briard-Bion: Writing -review & editing, Methodology; Fanny Guyomarc'h: Writing -review & editing, Methodology; Valérie Gagnaire: Writing -review & editing, Visualization, Methodology, Supervision.

Acknowledgments

N. Illikoud and F. Rosa do Carmo are recipients of post-doctoral fellowships from the French FUI (Fond Unique Interministériel) program, N °25, ProBioVeg. This study is part of the PhD project of F. Canon funded by the French National Research Institute for Agriculture, Food and Environment (INRAE) and the Brittany region. The authors offer their thanks to Victoria Chuat, from CIRM-BIA, for storing, checking and maintaining the bacterial strains.

was provided by the CIRM-BIA ( https://collection-cirmbia.fr/ ), international center for microbial resources, located in the STLO research unit ( https://www6.rennes.inrae.fr/stlo ). The work was performed at STLO, a research unit of: Institution: INRAE & Institut Agro City: Rennes

URL to data: https://www.ebi.ac.uk/pride/archive/projects/

Data Availability

Data from a proteomic comparative analysis highlight differential adaptation of Lactobacillus delbrueckii subsp. bulgaricus to cow milk versus to soy milk environments (Original data) (PRIDE).

Ethics Statements

This manuscript adheres to Ethics in publishing standards ( https://www.elsevier.com/ journals/data-in-brief/2352-3409/guide-for-authors ).

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Phuong Hoang Nguyen

Pierre Tufféry

Philippe Derreumaux

Dynamics of amyloid formation from simplified representation to atomistic simulations

Keywords: amyloid, aggregation, simulations, intrinsically disordered proteins, bulk solution, membranes, inhibitors, Aβ, tau, α-synuclein

Amyloid fibril formation is an intrinsic property of short peptides, non-disease proteins and proteins associated neurodegenerative diseases. Aggregates of the Aβ and tau proteins, the α-synuclein protein, and the prion protein are observed in the brain of Alzheimer's, Parkinson's and prion disease patients, respectively. Due to the transient short-range and long-range interactions of all species and their high aggregation propensities, the conformational ensemble of these devastating proteins, the exception being for the monomeric prion protein, remains elusive by standard structural biology methods in bulk solution and in lipid membranes. To overcome these limitations, an increasing number of simulations using different sampling methods and protein models have been performed. In this article, we first review our main contributions to the field of amyloid protein simulations aimed at understanding the early aggregation steps of short linear amyloid peptides, the conformational ensemble of the Aβ40/42 dimers in bulk solution, and the stability of Aβ aggregates in lipid membrane models. Then we focus on our studies on the interactions of amyloid peptides/inhibitors to prevent aggregation, and long amyloid sequences, including new results on a monomeric tau construct.

Introduction

Aβ, tau, α-synuclein and prion protein misfolding and aggregation in the central nervous system lead to three neurodegenerative diseases in elderly, Alzheimer's, Parkinson's and prion diseases (1-3). The senile extracellular Alzheimer's disease plaques between neurons contain Aβ42 and Aβ40 peptides, Aβ42 being more toxic than Aβ40, with Aβ sequence consisting of a charged hydrophilic N-terminus (residues 1-16), the central hydrophobic core (CHC, residues 17-21), a charged region (residues [START_REF] Lasagna-Reeves | Preparation and characterization of neurotoxic tau oligomers[END_REF][START_REF] Fusco | Structural basis of membrane disruption and cellular toxicity by αsynuclein oligomers[END_REF][START_REF] Graen | Transient Secondary and Tertiary Structure Formation Kinetics in the Intrinsically Disordered State of α-Synuclein from Atomistic Simulations[END_REF][START_REF] Nguyen | Structures of the intrinsically disordered Aβ, tau and αsynuclein proteins in aqueous solution from computer simulations[END_REF][START_REF] Santini | Pathway complexity of Alzheimer's betaamyloid Abeta16-22 peptide assembly[END_REF][START_REF] Santini | In silico assembly of Alzheimer's Abeta16-22 peptide into beta-sheets[END_REF][START_REF] Wei | Sampling the self-assembly pathways of KFFE hexamers[END_REF][START_REF] Melquiond | Following the aggregation of amyloid-forming peptides by computer simulations[END_REF] and a hydrophobic Cterminus (residues 30-40/42). The Aβ peptides result from the proteolytic cleavage of the transmembrane amyloid precursor protein (APP) by secretases (2). The senile intracellular Alzheimer's disease filaments in neurons are made of hyperphosphorylated tau proteins of 441 residues consisting of a charged N-terminal (residues 1-207), a proline-rich region, a microtubule-binding domain of four repeats spanning residues 244-368, and a C-terminal region (4). The α-synuclein protein has a membrane-binding domain (residues 1-66), a hydrophobic region spanning residues 67-96, and a highly charged C-terminal (residues 97-140) implicated in calcium binding and in protein-protein interactions at the surface of synaptic vesicles (1, 5). In contrast to the monomeric disordered Aβ, tau and α-synuclein structures, the monomeric prion structure obtained by nuclear magnetic resonance (NMR) using solution full-length and C-terminal domain from both mouse and Syrian hamster prion reveals an unstructured N-terminal fragment ) and a C-terminal globular domain (125-228) with three α-helices and a short antiparallel β-sheet between helices 1 and 2 (6). Despite low sequence identity, the four proteins aggregate over time to amyloid fibrils with a common cross-β structure (7-10).

Aggregation is a non-equilibrium process where the amyloid proteins form fibrils with sigmoidal kinetics profiles in which the proteins self-assemble (lag-phase) prior to fibril elongation (growth phase), followed by a plateau where the fibrils and free monomers are in equilibrium (saturation phase) (7, 11). The lag phase time of fibrillation and the fibril growth time vary with amino acid length, mutation, temperature, pH, concentration, agitation, shear forces, metal ions, crowding, and the presence of lipids such as cholesterol [START_REF] Meisl | Kinetic Analysis of Amyloid Formation[END_REF][START_REF] Buttstedt | Different morphology of amyloid fibrils originating from agitated and non-agitated conditions[END_REF][START_REF] Luo | Macromolecular crowding favors the fibrillization of β2-microglobulin by accelerating the nucleation step and inhibiting fibril disassembly[END_REF][START_REF] Xu | Intrinsic determinants of Aβ(12-24) pH-dependent self-assembly revealed by combined computational and experimental studies[END_REF][START_REF] Habchi | Cholesterol Catalyses Aβ42 Aggregation through a Heterogeneous Nucleation Pathway in the Presence of Lipid Membranes[END_REF]. Of note, in vitro aggregation of full-length tau requires the addition of heparin due to the highly charged N-and C-terminal regions and is impacted by the pattern of hyperphosphorylation (16).

The molecular mechanisms involved in amyloid fibril formation are well described by primary nucleation, (fragmentation and surface-catalyzed) secondary nucleation, and elongation. It is now possible to determine the rate constants of each of these microscopic processes operating at multiple time-scales by fitting the solutions of non-linear chemical master equations to the experimental aggregation curves (11, 15). Such a knowledge is essential to design molecular inhibitors that target selectively mature amyloid fibrils or intermediate oligomeric species and interfere with the desired microscopic aggregation step.

Interestingly, the Aβ, tau and α-synuclein proteins have hydrophobic aggregation-prone regions which by themselves form amyloid fibrils in vitro, e.g., the CHC consisting of 17 LVFFA 21 and the C-terminus 37 GGVIA 42 in Aβ, the PH6 and PH6* motifs consisting of PHF6*: 275 VQIINK 280 and PHF6: 306 VQIVYK 311 in tau, and the NAC region in α-synuclein (17-19). These fragments have therefore represented ideal models for understanding amyloid fibril formation by computational studies (20, 21).

It has also been demonstrated that the early formed aggregates of Aβ, tau and α-synuclein are toxic [START_REF] Hardy | The Amyloid Hypothesis of Alzheimer's Disease: Progress and Problems on the Road to Therapeutics[END_REF][START_REF] Lasagna-Reeves | Preparation and characterization of neurotoxic tau oligomers[END_REF][START_REF] Fusco | Structural basis of membrane disruption and cellular toxicity by αsynuclein oligomers[END_REF]. This finding has also motivated many computational studies to disclose the structures of the oligomeric species in various environments. This is a challenge computationally due to the size and heterogeneity of the conformational ensemble to be explored and the need of accurate force fields to model intrinsically disordered proteins (24, 25). The next sections review our main contributions to the field of amyloid protein simulations. We also present new results of the application of the PEP-FOLD framework on a monomer of tau construct.

Early aggregation steps of short linear amyloid peptides in bulk solution

To accelerate aggregation, all simulations use a peptide concentration of mM, namely three and six orders of magnitude higher than in vitro and in vivo conditions, respectively. Various oligomers of the Aβ16-22, Aβ37-42, Aβ25-35, Aβ11-25, KFFE, GNNQQNQ and NNQQ from yeast protein Sup35, NFGAIL and SNNFGAILSS from islet amyloid polypeptide (IAPP), and NHVTLSQ from β2-microglobulin peptides (Table 1) were explored by coarse-grain (OPEP) and atomistic simulations, based on molecular dynamics (MD), Monte Carlo (MC), the Activation-Relaxation Technique (ART), standard replica exchange MD (REMD) and replica exchange with solute tempering (REST2) . OPEP, which consists of an all-atom representation of the backbone, one bead per side-chain, the exception being proline, and an implicit solvent model, was parametrized on well-folded peptides and proteins (47-52). ART accelerates the search of lowest-energy microstates by finding repeatedly activated mechanisms that connect minima via first-order saddle-points and accepting/rejecting them by the Metropolis criterion (53-55).

Independently of the protein (coarse-grain or atomistic) and solvent (implicit or explicit) representation and the amino acid sequence, self-assembly of up to 20 peptides starts by a hydrophobic collapse and the formation of disordered oligomers, which involve in time to transient β-rich aggregates (35-37). These heterogeneous β-rich assemblies can have various sheet-to-sheet pairing angles ranging from in parallel to perpendicular, and form open and closed β-barrels (28-30, 37, 39, 42). Of note, these topologies have been observed in vitro by X-ray crystallography of macrocyclic β-sheet mimics (56) and designed amyloid peptides (57). For all peptides, however, the free energy landscape is dominated by amorphous aggregates [START_REF] Mo | Structural diversity of the soluble trimers of the human amylin(20-29) peptide revealed by molecular dynamics simulations[END_REF][START_REF] Nasica-Labouze | A multiscale approach to characterize the early aggregation steps of the amyloid-forming peptide GNNQQNY from the yeast prion sup-35[END_REF][START_REF] Lu | Structural, thermodynamical, and dynamical properties of oligomers formed by the amyloid NNQQ peptide: insights from coarse-grained simulations[END_REF]. Our studies, as other computational studies [START_REF] Matthes | An Atomistic View of Amyloidogenic Self-assembly: Structure and Dynamics of Heterogeneous Conformational States in the Pre-nucleation Phase[END_REF][START_REF] Levine | Simulations of disordered proteins and systems with conformational heterogeneity[END_REF] highlight the impact of the monomeric amyloid-competent state, which corresponds to an extended conformation for short amyloid sequences, on (i) the kinetics of aggregation, e.g., the association/dissociation times of all oligomers and the conversion time between architectures, and (ii) on the thermodynamics of aggregation, such as the population of the substates having different β-strand mismatch (mixed parallel/antiparallel β-strands), βsheet sizes and topologies [START_REF] Santini | In silico assembly of Alzheimer's Abeta16-22 peptide into beta-sheets[END_REF][START_REF] Lu | Thermodynamics and dynamics of amyloid peptide oligomerization are sequence dependent[END_REF]. The simulations also emphasize a complex free energy landscape with βrich oligomers having in-register and out-of-register conformations and the existence of many kinetic traps

It is not yet understood which recent modern force field is the most appropriate for describing the aggregation pathways and the equilibrium ensemble of intrinsically disordered proteins (33, 62-66).

Similarly, the nucleus size (N*) for primary nucleation, corresponding to the highest free energy aggregate from which fibrils can grow, remains elusive, simulations predicting N* between 5 and 40 [START_REF] Nasica-Labouze | Amyloid beta Protein and Alzheimer's Disease: When Computer Simulations Complement Experimental Studies[END_REF][START_REF] Baftizadeh | Nucleation process of a fibril precursor in the C-terminal segment of amyloid-β[END_REF][START_REF] Šarić | Kinetics of spontaneous filament nucleation via oligomers: Insights from theory and simulation[END_REF][START_REF] Lee | Mechanisms and rates of nucleation of amyloid fibrils[END_REF][START_REF] Nguyen | Understanding amyloid fibril nucleation and aβ oligomer/drug interactions from computer simulations[END_REF][START_REF] Tran | Lattice model for amyloid peptides: OPEP force field parametrization and applications to the nucleus size of Alzheimer's peptides[END_REF][START_REF] Chiricotto | Coarsegrained and all-atom simulations towards the early and late steps of amyloid fibril formation[END_REF][START_REF] Sterpone | Multiscale simulations of biological systems using the OPEP coarse-grained model[END_REF]). Yet, this knowledge is important as nucleation seed size determines amyloid clearance and establishes a barrier to prion appearance in yeast (75). Another important issue, when performing coarse-grained simulations with an implicit solvent representation, is that exchange of momenta between the peptide's particles and the solvent and solvent mediated correlations are ignored. In this context, we coupled Lattice-Boltzmann MD simulation, which includes naturally hydrodynamics interactions, to the CG OPEP model [START_REF] Sterpone | Protein Simulations in Fluids: Coupling the OPEP Coarse-Grained Force Field with Hydrodynamics[END_REF][START_REF] Chiricotto | Multiscale simulation of molecular processes in cellular environments[END_REF] and performed simulations of Aβ16-22 peptides in bulk solution (78-80). Two main results are important. For a system of 100 peptides, hydrodynamics interactions augment the oligomer size of the first two clusters and the exchange of peptides compared to the results of standard Langevin dynamics (78). For a system of 1000 peptides at a concentration of 60mM, we can follow the formation and growth of a large elongated oligomer and its slow β-sheet structuring resulting from the fusion and dissociation of small disordered aggregates. Many mechanisms are observed from elongation to surface-catalyzed effect, and in particular the lateral growth mechanism on the surface of prefibrillar states, which is sustained by long-range hydrodynamics correlations, and allows the formation of large branched structures consisting of 600 peptides, spanning a few tens of nanometers and hosting annular pores of dimensions 3-5 nm (79, 80). This computational simulation at a quasi-atomistic representation of a system of unprecedented size, previous simulations exploring at most 125 peptides with very few degrees of freedom, illustrates the critical contribution of secondary nucleation to amyloid aggregation kinetics. It is interesting that amyloid plaques with annular pores have been evidenced by electron microscopy and antibodies in the brain of Alzheimer's disease patients (81). Of note our early phase of aggregation would approach the second time scale at nM (in vivo) concentrations.

Dimeric Aβ40/42 states in bulk solution

The wild-type (WT) Aβ40/42 dimers, known to be structurally heterogeneous by experiments, are of high interest as they are the smallest species to lead to tau hyperphosphorylation and neuritic degeneration, and their levels have been found to increase sharply and correlate with plaque load in brain tissue of Alzheimer's disease patients (82). Experiments have shown that the English (H6R) familial Alzheimer disease mutation, the Tottori (D7N) mutation, the Flemish (A21G) mutation and the Iowa (D23N) mutation speed up the fibril formation process of both Aβ40 and Aβ42 peptides in vitro and increase toxicity to cells (21). In contrast, the engineered disulfide-bond-locked double (S8C) mutant has been shown to form an exclusive homogeneous and neurotoxic dimer (83). Mutations at position 2 have dramatic impact on AD risk; A2V is causative, A2T is protective. The A2V mutation enhances aggregation kinetics while the A2T mutation only retards amyloid fibril by increasing the lag phase time. Interestingly the mixture of WT and A2V also retards fibril formation and protects against Alzheimer's disease (84-86).

The Aβ40 and Aβ42 alloforms were subject to atomistic REMD in explicit solvent (87, 88). Using several force fields, the equilibrium configurations are found disordered, with cross-collision sections, hydrodynamics radii, and SAXS profiles of the ensembles independent of the force field. Intramolecular βhairpins spanning residues 17-21 and 30-36 are however observed with a population varying from 1.5 to 13% according to the force field. The ensemble of both alloforms is stabilized by nonspecific interactions with many hydrophobic residues exposed to solvent, explaining therefore their propensity to be toxic at the dimeric level. Simulations also reveal that the Aβ42 dimer has a higher propensity than the Aβ40 dimer to form β-strands at the CHC and at the C-terminal (residues 30-40), consistent with other computational studies (89,90). The dimers have no defined interfaces, and the random organization with transient secondary structures is substantially preferred over two chains with β-hairpin and fibril-like conformations (88). The formation and the transition between the two latter β-rich dimers have also been discussed by other studies (91, 92).

Using atomistic simulations, the dimers of S8C and WT Aβ42 have the same secondary structures and cross-collision sections. Upon S8C mutation, the lifetime of the intramolecular three-stranded spanning residues 17-21, 30-36, and 39-41 is increased by a factor of 3. This single common structural feature shared by both species, which does not exist in the Aβ40 WT dimers, is likely to contribute to Aβ42 toxicity (93). The H6R, D7N, A21G, D23N mutations change the population of β-strand at the CHC region and at the Cterminus, and impact the population of intramolecular and intermolecular salt bridges involving E22, D23 and K28 by reducing the formation time of the loop region for a fibril-like conformation. We also find, consistent with experimental studies, that these mutations produce different effects on Aβ dimer depending on whether they occur in Aβ40 or Aβ42 (90, 94-96).

Comparing atomistic REMD simulations of WT, WT-A2V and A2V-A2V Aβ40 dimers, we find that upon single mutation, the intrinsic disorder and the intermolecular potential energies are reduced and the population of intramolecular three-stranded β-sheets is increased (97). A reduced intrinsic disorder upon A2V mutation was already discussed for the Aβ28 monomer (98). Analyzing REMD simulations of WT-A2T Aβ40 dimer (99), we provide evidence that the retard in the lag phase time of fibrillation results from an increase of intrapeptide stability over interpeptide stability in the heterozygous dimers. This finding was further corroborated by other simulations of A2T and A2V Aβ42 dimers (100, 101). Investigating the local structure and dynamics of hydration, it was also found that the survival probability of order water molecules decays more rapidly for the Aβ N-terminus AD causative (A2V, H6R, D7N and D7H) mutants than the A2T protective mutant (102).

4.Oligomeric states of Aβ in lipid membrane models

Interactions of amyloid oligomers with cell membranes are believed to contribute to toxicity. On one hand, the membrane can trigger amyloid fibril formation at a lower peptide concentration (103, 104). On the other hand, accumulation of oligomers on the membrane surface can impart inequal stress on bilayer, extract lipids into and contribute to the formation of stable phospholipids/oligomers complexes or create pores transporting Ca 2+ ions through the membrane (105-110). Due to structural heterogeneity of the oligomers, many pores of various inner diameters made of different number of Aβ subunits have been modelled based on atomic force microscopy images and MD simulations (110).

Based on size exclusion chromatography, transmission electron microscopy, circular dichroism and NMR information on the structural environment experienced by the Q15, N27 and M35 side chains (111), we recently modeled a tetrameric β-barrel consisting of two distinct β-hairpins, with an assymetric arrangement of eight antiparallel β-strands and an inner pore diameter pore of 0.7 nm for both alloforms of Aβ (112).

Using extensive atomistic REMD simulations, we found that this barrel exists transiently for Aβ42 and not for Aβ40 within a dipalmitoylphosphatidylcholine (DDPC) lipid bilayer membrane and this may explain the higher toxicity effect of Aβ42 than its Aβ40 counterpart (112). The same simulations indicate that the lower and higher induced toxicity of the A2T and D23N mutants cannot be correlated to their tetrameric β-barrel pore-forming probabilities at least in a DDPC membrane environment (113), but addition of cholesterol in the membrane composition could make a difference (114). As the structures of amyloid oligomers inserted into a membrane remain elusive, we explored the stability of Aβ11-42 trimers with parallel (U-shape fibril) and antiparallel (β-hairpin) β-sheet structures in a DPPC membrane. Our REMD simulations strongly suggest that these two assemblies represent minimal seeds or nuclei for the formation of amyloid fibrils, a variety of β-barrel pores and various aggregates for Aβ sequences (115). Apart from β-barrel pores, α-helical pores are also possible (107, 116). The equilibrium structures are in all cases dependent, however, on the lipid composition (117). For instance, it is known that the omega-3 polyunsaturated fatty acid slows the progression of Alzheimer's disease, while its omega-6 counterpart is linked to increased risk. We showed by MD simulations that variation in the abundance of the 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), omega-3 and omega-6 modulates the conformational ensemble of the Aβ29-42 dimer (118).

Study of inhibitors/Aβ amyloid oligomers

Understanding the mechanistic determinants of Aβ amyloid-inhibitor interaction is continuously pursued to develop more efficient drugs (70). A plethora of inhibitors aimed at targeting either oligomers, the secondary nucleation or fibril elongation have been tested in bulk solution and in vivo (119). We performed numerous simulations to study the detailed interactions of Aβ fragments, Aβ40 and Aβ42 oligomers with a large set of inhibitors including N-methylated peptide, EGCG, curcumin, resveratrol, 1,4-naphthoquinon-2-yl-L-tryptophan (NQ-Trp), astaxanthin and betanin (120-127). We have shown, as other computational studies (128-130), that there are many binding sites with small occupancies and contact surfaces. Even with a binding pocket, there are multiple binding modes, demonstrating the transient character of Aβ oligomers/inhibitors interactions in bulk solution. This is consistent with the absence of non-specific interactions as evidenced by NMR and the low affinity of drugs for Aβ monomers and small oligomers (131) .

Besides small molecules, carbon nanoparticles such as fullerene and carbon nanotubes can also impede the fibrillation of Aβ and β2m proteins (132, 133). In this context, we have shown by atomistic REMD simulations that carbon nanotubes impact both the primary and secondary nucleation mechanism by destabilizing β-rich oligomers and fibril assemblies for both the Aβ16-22 and Aβ25-35 peptides (134, 135).

Several reasons have been put forward to explain the repetitive failures of small drugs and antibodies targeting Aβ oligomers (131).

As the lipid membrane can catalyze the formation of toxic amyloid intermediates, the interplay between amyloid oligomers, inhibitors and the lipid membrane should also be considered more systematically (136, 137). Similarly, clinical trials might consider a synergy of multiple drugs targeting Aβ and tau along with the use of laser or bubble cavitation techniques to destabilize amyloid aggregates (138-141).

Long amyloid constructs

The conformational ensemble of several monomeric prion fragments has been explored by coarsegrained OPEP Monte Carlo (48), and atomistic MD and REMD simulations. The CG simulations on PrP(143-158), starting from random states, showed that this fragment forms helix by itself consistent with experiments (142) . In contrast, the fragment PrP(128-164) was found to code either for an alpha/β topology, as found in the NMR structure of recombinant full-length PrP, or a β-hairpin spanning residues 142-167, as found by NMR for PrP(142-167) (142). Using atomistic MD, we showed that the helix H1 (residues 143-158) is rather stable upon P102L, M129V and G131V mutation and deletion of the residues coding the first or the second β-strand (143, 144); using REMD simulations, we found an intermediate state characterized by a significant detachment of helix H1 from PrP-core (145), consistent with other studies (146), which forms very easily for the E211D mutation by standard MD simulations (147). This is of interest as mutation at position 211 drives a switch between Creutzfeldt-Jakob disease (CJD) and Gerstmann-Sträussler Scheinker syndrome (147). Apart from the implication of H1 detachment in the early steps of aggregation, we found that the CJD-causing T183A mutation accelerates the conversion of the helix H2 into β-sheets in the monomer and dimer of PrP (148).

In comparison to the Aβ protein, the number of computational studies on alpha-synuclein monomer and long tau constructs is very small. We recently performed atomistic REMD simulations on the tau R3-R4 domain (residues 306-378) starting from the fibril topology (149). Note that the cryo-electron microscopy of full-length tau fibril in the brain of an individual with Alzheimer's disease reveals a R3-R4 domain with a Cshaped topology of eight β-sheets, while the rest of the protein is flexible (9). We found that the WT R3-R4 dimer populates elongated, U-shaped, V-shaped and globular topologies rather than C-shaped forms. MD simulations revealed that upon phosphorylation of Ser356 (pSer356) there is a substantial decrease of intermediates near the fibril like conformers compared to its WT counterpart (149). This result explains why WT K18 develops seeding activity more rapidly than pSer356 K18 (150).

Recently, we applied the chain-growth PEP-FOLD approach, developed for ab initio structure prediction of well-folded proteins (151-155) and protein-peptide complexes (156) to the monomeric αsynuclein protein (25). As expected from experiments, the α-synuclein monomer was found highly dynamic (Rg distribution varying 1.5 nm and 2.5 nm) and without any well-defined structure (65% of turn/coil). We observed however a high propensity of broken helices in the N-terminus consistent with α-synuclein's membrane binding properties (25) . monomer, the R3-R4 domain being known to be the core of the amyloid fibril. As seen on the free energy landscape, there are no dominant structures, the Rg distribution varying between 1.2 and 2.0 nm and the Cα end-to-end distance between 0.5 and 5.0 nm (Figure 1). Yet, β-strands signals are detected by the structural alphabet profile along the sequence, and regions with a β-strand propensity of more than 0.2 encompass residues 306-315, 317-321, 329-332, 335-345, 349-354, 361-364 and 367-378 (Figure 2). It is interesting that these regions fit well the β-strands observed in the fibril that encompass residues 306-311, 313-322, 327-331, 336-341, 343-347, 349-354, 356-363 and 368-378 (9). Most of structures are however devoid of βsheets and are random coil with some helical content (Figures 1A-1C), and a few structures have the propensity to form a double U-shape free of intramolecular H-bonds (Figure 1D). Overall, we do not find any evidence of a fibril-like conformation encoded in the dimeric ensemble.

Conclusions

We have reviewed our contributions to the field of amyloid simulations. Conformational ensembles of amyloid monomer and oligomers have been addressed from bulk solution, membrane environment to complexes with different classes of inhibitors. Improved coarse-grained and atomistic force fields should tell us more about the link between oligomers and toxicity. This requires however to mimic in vivo conditions such as crowding, and lipids, and to play with aging parameters such as variation of the fluid flow in the brain extracellular space (157, 158).

(152) Maupetit, J., Derreumaux, P. and Tufféry, P. (2010) A fast method for large-scale de novo peptide and miniprotein structure prediction. J Comput Chem., 31 : 726-738. ---------------------------------------------------------------------------------------------------------------------------------------------------1 Simulations in bulk solution unless specified. 2 Aβ16-22 dimer and trimer, Aβ11-25 tetramer, NFGAIL dodecamer, KFFE dimer up to heptamer, and 3-mers, 12-mers and 20-mers of GNNQQNY. 3 All-atom for Aβ16-22 dimer, OPEP for 4-mers up 16-mers of β2m [START_REF] Müller-Schiffmann | Molecular Engineering of a Secreted, Highly Homogeneous, and Neurotoxic Aβ Dimer[END_REF][START_REF] Rousseau | The Alzheimer disease protective mutation A2T modulates kinetic and thermodynamic properties of amyloid-β (Aβ) aggregation[END_REF][START_REF] Zheng | Amyloid β-Protein Assembly: Differential Effects of the Protective A2T Mutation and Recessive A2V Familial Alzheimer's Disease Mutation[END_REF][START_REF] Messa | The peculiar role of the A2V mutation in amyloid-β (Aβ) 1-42 molecular assembly[END_REF][START_REF] Tarus | Structures of the Alzheimer's Wild-Type Aβ1-40 Dimer from Atomistic Simulations[END_REF][START_REF] Man | High-Resolution Structures of the Amyloid-β 1-42 Dimers from the Comparison of Four Atomistic Force Fields[END_REF][START_REF] Côté | Distinct Morphologies for Amyloid Beta Protein Monomer: Aβ1-40, Aβ1-42, and Aβ1-40(D23N)[END_REF], Aβ [START_REF] Sibille | Structural impact of heparin binding to full-length Tau as studied by NMR spectroscopy[END_REF][START_REF] Balbach | Amyloid fibril formation by A beta 16-22, a seven-residue fragment of the Alzheimer's beta-amyloid peptide, and structural characterization by solid state NMR[END_REF][START_REF] Inouye | Structure of core domain of fibril-forming PHF/Tau fragments[END_REF][START_REF] Bodles | Identification of the region of non-Abeta component (NAC) of Alzheimer's disease amyloid responsible for its aggregation and toxicity[END_REF][START_REF] Mousseau | Exploring the early steps of amyloid peptide aggregation by computers[END_REF][START_REF] Nasica-Labouze | Amyloid beta Protein and Alzheimer's Disease: When Computer Simulations Complement Experimental Studies[END_REF][START_REF] Lasagna-Reeves | Preparation and characterization of neurotoxic tau oligomers[END_REF] 8-mers and IAPP(20-29) trimer. 4 All-atom for Aβ16-22 dimer and trimer, 7-mers of β2m(83-89) and 16-mers of Aβ37-42; OPEP for ccβ peptide, 2-mers and 3-mers of Aβ16-22, 6-mers of Aβ25-35, 20mers of NNQQ and 3-mers. 5 12-mers and 20-mers of GNNQQNY.

( 26 -

 26 31, 40, 42), as it has been observed experimentally in the late steps of aggregation (60) or computationally by Markov chain models in the lock phase of oligomers (61).

( 153 )

 153 Thévenet, P., Shen, Y., Maupetit, J., Guyon, F.,Derreumaux, P. and Tufféry, P. (2012) PEP-FOLD: an updated de novo structure prediction server for both linear and disulfide bonded cyclic peptides. Nucleic Acids Res., 40 : W288-W293. (154) Shen, Y., Maupetit, J., Derreumaux, P. and Tufféry, P. (2014) Improved PEP-FOLD Approach for Peptide and Miniprotein Structure Prediction. J Chem Theory Comput., 10 : 4745-4758. (155) Sutherland, G.A., Grayson, K.J., Adams, N.B.P. et al. (2018) Probing the quality control mechanism of the Escherichia coli twin-arginine translocase with folding variants of a de novo-designed heme protein. J Biol Chem., 293 : 6672-6681. (156) Lamiable, A., Thévenet, P., Rey, J.,, Vavrusa, M., Derreumaux P. and Tufféry, P. (2016) PEP-FOLD3: faster de novo structure prediction for linear peptides in solution and in complex. Nucleic Acids Res., 44 : W449-W454. (157) Ngo, S.T., Nguyen, P.H. and Derreumaux, P. (2021) Cholesterol molecules alter the energy landscape of small Aβ 1-42 oligomers. J Phys Chem B., doi : 10.1021/acs.jpcb.1c00036 (158) Ramamoorthy, A.; Sahoo, B. R.; Zheng, J.; Chiricotto, M.; Straub, J. E.; Dominguez, L.; Shea, J-E.; Dokholyan, N.V.; De Simone, A. et al. (2021) Amyloid Oligomers: A Joint Experimental/Computational Perspective on Alzheimer's Disease, Parkinson's Disease, Type II Diabetes and Amyotrophic Lateral Sclerosis. Chem. Rev., Chem Rev., doi: 10.1021/acs.chemrev.0c01122.

Figure 1 .

 1 Figure 1. PEP-FOLD free energy landscape (in Kcal/mol) of the tau (306-378) monomer projected on the radius of gyration and the end-to-end distance. The N-terminus is denoted by a sphere, and the C-terminus by a square.

Figure 2 :

 2 Figure 2: Structural alphabet profile predicted for tau-306-378. Green tones correspond to extended conformations, red tones to helical conformations and blue one to coil regions. Extended regions. Each column corresponds to a fragment of 4 amino-acids (first column corresponds to VQIV, last column corresponds to KLTF).

Table 1 .

 1 Summary of the different amyloid systems studied.

	Applications 1	Sampling	Force Field	Ref.
	Aggregation of short linear peptides	ART 2	OPEP	26-32
		MD 3	OPEP, all-atom 33-37
		REMD 4	all-atom	38-40
			OPEP	41-46
		REST2 5	all-atom	42
	Nucleation sizes of Aβ16-22 and Aβ37-42	REMC	on-lattice OPEP 71-73
	Impact of hydrodynamics on early
	aggregation steps	LBMD	OPEP	78,79
	Aβ40 WT, Aβ42 WT dimers	REMD	all-atom	87,88
	Aβ40 WT, Aβ42 WT, D23N monomers and dimers ART, REMD	OPEP	89,90
	Aβ16-35 monomer and dimer	REMD	OPEP	92
	Aβ42 (S8C) dimers	REMD	all-atom	93
	Aβ42 H6R, D7N dimers	MD	all-atom	94-96
	Aβ42 A2T, A2V dimers	REMD	all-atom	97,99
	Aβ1-28 WT and A2V monomers	REMD	all-atom	98
	Trimeric Aβ11-42 structures in DPPC membrane	REMD	all-atom	115
	Tetrameric Aβ40, Aβ42 WT, D23N, A2T β-barrels
	and Aβ helix bundles in DPPC lipid membrane	REMD	all-atom	112,113,116
	Aβ29-42 dimer with membranes containing
	omega-3 and omega-6 fatty acids	MD	all-atom	118
	Aβ1-28 / NQ-Trp	REMD	all-atom	124,125
	Aβ16-22 oligomers/N-methylated inhibitor	MD	OPEP	120
	Aβ17-42 trimers / small molecule inhibitors	REMD	all-atom	121
	Aβ40/42 dimers / small molecule inhibitors	MD, REMD	all-atom	122,123,127
	Aβ16-22, Aβ25-35 oligomers / carbon nanotubes MD, REMD	all-atom	134,135
	prion (127-164) monomer	MC	OPEP	142
	prion WT monomer	MD	all-atom	143,144,147
		REMD	all-atom	145
	prion (125-228) dimer	MD	OPEP	148
	Dimeric tau (306-378) construct
	and phosphorylation	REMD, MD	all-atom	149
	α-synuclein monomer	PEP-FOLD	OPEP	25
	Monomeric tau (306-378)	PEP-FOLD	this work
	-

ACKNOWLEDGMENTS

We acknowledge support by the "Initiative d'Excellence" program from the French State (Grant "DYNAMO", ANR-11-LABX-0011-01, and "CACSICE", ANR-11-EQPX-0008).

AUTHOR INFORMATION

Corresponding Author E-mail: philippe.derreumaux@ibpc.fr ORCID Phuong Hoang Nguyen: 0000-0003-1284-967X Pierre Tufféry : 0000-0003-1033-9895 Philippe Derreumaux: 0000-0001-9110-5585

Notes

The authors declare no competing financial interest.