Carla Palacios-Gorba

Alexandra Moura

Yuval Markovich

Nathalie Tessaud-Rita

Ángel Gómez-Martín

Hélène Bracq-Dieye

Jesús Gomis

Guillaume Vales

María Pastor-Martín

Pierre Thouvenot

Carles Escrig

Alexandre Leclercq

Marc Lecuit

Juan J Quereda
email: juan.quereda@uchceu.es.

Genomic characterization of Listeria spp. isolated from tonsils, udder and feces of domestic dairy ruminants in Spain

Keywords: Listeria spp., dairy, ruminant, tonsil, udder, WGS

à 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

Listeriosis is recognized as a significant public health issue due to the severity of the disease and its mortality rate (20% to 30%) [START_REF] Charlier | Clinical features and prognostic factors of listeriosis: the MONALISA national prospective cohort study[END_REF][START_REF] Quereda | Pathogenicity and virulence of Listeria monocytogenes: A trip from environmental to medical microbiology[END_REF]. Only two species of the genus Listeria, Listeria monocytogenes (Lm) and Listeria ivanovii (Liv), are considered facultative intracellular pathogens [START_REF] Quereda | Pathogenicity and virulence of Listeria monocytogenes: A trip from environmental to medical microbiology[END_REF]. Lm can cause serious invasive disease (including central nervous system and fetal-placental infection) in animals and humans [START_REF] Charlier | Clinical features and prognostic factors of listeriosis: the MONALISA national prospective cohort study[END_REF][START_REF] Quereda | Pathogenicity and virulence of Listeria monocytogenes: A trip from environmental to medical microbiology[END_REF]. Although Liv is mainly an animal pathogen, its isolation from human patients with listeriosis symptoms has been reported [START_REF] Guillet | Human listeriosis caused by Listeria ivanovii[END_REF]. Listeria innocua (Lin) is a close relative of Lm and is nonpathogenic to mammals, although atypical virulent isolates have been reported [START_REF] Moura | Atypical Hemolytic Listeria innocua Isolates Are Virulent, albeit Less than Listeria monocytogenes[END_REF].

Healthy people and domestic ruminants tend to be asymptomatic carriers that shed Lm, Liv and non-pathogenic Listeria spp. in their feces [START_REF] Hafner | Listeria monocytogenes faecal carriage is common and depends on the gut microbiota[END_REF][START_REF] Palacios-Gorba | Ruminant-associated Listeria monocytogenes isolates belong preferentially to dairy-associated hypervirulent clones: a longitudinal study in 19 farms[END_REF]. Fecal shedding of Lm could lead to contamination of agricultural environments and raw produce which are then transported to food processing industries [START_REF] Walland | Listeria monocytogenes infection in ruminants: Is there a link to the environment, food and human health? A review[END_REF]. In dairy farms, Lm is usually transmitted to bulk tank milk as an outcome of environmental or fecal contamination of the udder [START_REF] Castro | Occurrence, Persistence, and Contamination Routes of Listeria monocytogenes Genotypes on Three Finnish Dairy Cattle Farms: a Longitudinal Study[END_REF].

While the incidence of listeriosis in humans is low compared with other foodborne infections, listeriosis is recognized as a significant public health issue due to its severity and high mortality rate [START_REF] Charlier | Clinical features and prognostic factors of listeriosis: the MONALISA national prospective cohort study[END_REF]. Lm population is diverse, comprising 4 lineages (ref), and its population structure can be determined by using multilocus sequence typing (MLST) and core genome MLST (cgMLST). MLST, which analyzes seven housekeeping genes, groups Lm isolates in sequence types (STs) and clonal complexes (CCs or clones) (Ref Ragon2008). cgMLST, which analyzes 1748 core loci, allows Lm classification in cgMLST types (CTs) and sublineages (SLs) [START_REF] Moura | Whole genome-based population biology and epidemiological surveillance of Listeria monocytogenes[END_REF]. Listeria spp. surveillance in animal production settings normally targets fecal or environmental samples, however, information on Listeria spp. presence in organs such as tonsils and udders is still scarce. Moreover, there is no information about asymptomatic carriage of Listeria spp. in ruminants' organs without fecal shedding. Here, we used whole-genome sequencing to subtype and characterize Listeria spp. isolated from feces, tonsils, and udders from domestic dairy ruminants. Increasing the understanding of the reservoirs and genomic characteristics of Lm and Liv will help to improve their control, enhancing human and animal disease prevention.

Materials and methods

Animals and sampling

In the time frame of December 2019 to May 2021, 93 dairy sheep (Ovis orientalis aries breed Lacaune, Manchega, Assaf, and crossbreed), 10 dairy goats (Capra aegagrus hircus breed Guadarrama and crossbreed), and 183 dairy cattle (Bos primigenius taurus breed Holstein-Friesian) were slaughtered in two slaughterhouses (slaughterhouse A located in Madrid for small ruminants and slaughterhouse B located in Valencia for cattle). Each slaughtering was followed by a systematic cleaning and disinfection protocol. Veterinary supervision was provided throughout the slaughter and Good Manufacturing Practices (GMP), Sanitation Standard Operating Procedures (SSOP), Hazard Analysis Critical Control Point (HACCP), and the traceability was strictly followed. Additionally, 29 dairy sheep and 1 dairy goat were culled for educational purposes to the authors' university. An intravenous injection of barbiturate was used to euthanize these additional animals, and complete necropsies were carried out immediately following death. Sheep and goats originated from 36 farms located in the provinces (administrative division in Spain) from Valencia, Albacete, Cuenca, Toledo, Ciudad-Real, and Madrid (Table S1). Cattle originated from 13 farms located in the provinces of Murcia, Albacete, Teruel, Alicante and Valencia (Table S1 and Fig. S2). Approximately 10 g of tonsils and udders were removed aseptically and stored in sterile bags. Tonsils were sampled while the carcasses were hanging upside down which implies that surface contamination was not possible. Udders were sampled once they were removed from the animal and located on a tray. A sample from the inner part of the udder was obtained by applying a cut on a previous decontaminated part of the organ surface with alcohol. Veterinary inspectors assumed that slaughtered animals were healthy based on the lack of abnormal behavior, clinical signs and macroscopic abnormalities in the thoracic and abdominal cavities of the carcasses. Thoracic and abdominal organs were inspected to evaluate position, size, and color.

Serosal surfaces were examined, and organs evaluated for any grossly visible abnormalities in structure and/or pathologic processes. The mandibular, the lung (cervical, mediastinal, tracheobronchial) and mesenteric lymph nodes, were also identified and evaluated. Collected feces from the rectum of each ruminant were stored in independent sterile bags. Each animal sample (tonsils, udder, and rectum) was handled with a new and independent sterile scalpel to prevent contamination. All samples were collected using disposable gloves under aseptic conditions. Within 24h, the entire sample material was transported to the laboratory at 5+/-3°C.

Listeria species isolation and identification.

Listeria spp. were isolated as previously described [START_REF] Palacios-Gorba | Ruminant-associated Listeria monocytogenes isolates belong preferentially to dairy-associated hypervirulent clones: a longitudinal study in 19 farms[END_REF][START_REF] Quereda | Listeria valentina sp. nov., isolated from a water trough and the faeces of healthy sheep[END_REF][START_REF] Palacios-Gorba | Listeria spp. Isolated from Tonsils of Wild Deer and Boars: Genomic Characterization[END_REF]. Briefly, 8 g of rectal fecal samples, 5 g of tonsils, and udders were disrupted and diluted 1/10 in Half-Fraser broth (Scharlab, Spain). An enrichment step of 24 h at 30 °C was performed on homogenized samples. Incubated suspensions were transferred into Fraser broth (Scharlab, Spain), and incubated for 24 h at 37 °C, before being transferred to RAPID'L. mono plates (BioRad, USA) and incubated at 37 °C for 24 h. RAPID'L. mono plates support growth of L. monocytogenes (PIPLC+/xylose-, blue colonies), L. ivanovii (PIPLC+/xylose+,blue-green colonies with yellow halo), and other Listeria spp. (white colonies). When more than one type of colony was present in RAPID'L.mono plates, one blue (no plate contained simultaneously L. monocytogenes and L. ivanovii) and one white Listeria spp.

were picked and further examined in selective Oxford agar plates for Listeria (Scharlab, Spain). From each positive sample, one isolate colony was obtained (if in the same plate both pathogenic and non-pathogenic Listeria spp. were present, one colony of each was selected) and sent to the World Health Organization Collaborating Center Listeria (Institut Pasteur, Paris) for characterization. As described previously, Listeria isolates were identified using matrix-assisted laser desorption ionization-time of flight mass spectrometry applying the MicroFlex LT system with the last MBT library DB-7854 (Bruker Daltonics, Germany) [START_REF] Thouvenot | MALDI-TOF mass spectrometry-based identification of Listeria species in surveillance: A prospective study[END_REF] and by whole genome sequencing [START_REF] Quereda | Listeria valentina sp. nov., isolated from a water trough and the faeces of healthy sheep[END_REF].

DNA extraction, genome sequencing and sequence analysis.

DNA was purified with the NucleoSpin Tissue DNA kit (Machenry-Nagel) from 900 µL bacterial cultures grown overnight in 5 mL of brain heart infusion broth (Becton Dickinson, USA) at 37 °C.

DNA libraries were prepared with Nextera XT DNA Sample kit (Illumina, USA), following manufacturer's protocol, and run in NextSeq 500 platform (Illumina, USA) to obtain paired-end reads (2 x 150 bp run). Raw reads were trimmed with fqCleaner v.21.10 and assembled with SPAdes v.3.15.0, using automatic k-mers selection and the --only-assembler and -careful parameters.

In silico typing was performed from the assemblies using BIGSdb-Listeria v.1.33.1 [START_REF] Jolley | BIGSdb: Scalable analysis of bacterial genome variation at the population level[END_REF] using the genoserogrouping, MLST, cgMLST profiles, resistance and virulence schemes implemented within [START_REF] Moura | Whole genome-based population biology and epidemiological surveillance of Listeria monocytogenes[END_REF], as previously described [START_REF] Palacios-Gorba | Ruminant-associated Listeria monocytogenes isolates belong preferentially to dairy-associated hypervirulent clones: a longitudinal study in 19 farms[END_REF][START_REF] Moura | Whole genome-based population biology and epidemiological surveillance of Listeria monocytogenes[END_REF]. A maximum phylogeny was inferred from the core genome alignment built using Parsnp, v1.5.6 (as implemented in Harvest suite 1.3 [START_REF] Treangen | The harvest suite for rapid core-genome alignment and visualization of thousands of intraspecific microbial genomes[END_REF]), using IQ-TREE/2.0.6 [START_REF] Nguyen | IQ-TREE: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies[END_REF]. Trees were visualized with iTol v.4.2 [START_REF] Letunic | Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees[END_REF].

Statistical analysis

IBM SPSS Statistics version 25 was used for all statistical analyses . Statistical tests were considered significant at P<0.05. Chi-square (ꭕ²) tests were performed to determine the effect of geographical location of farms and animal species on the number of Lm or Liv carriage.

Data availability.

All sequences are publicly available at the European Nucleotide Archive (BioProject PRJEB51772) and BIGSdb-Listeria (https://bigsdb.pasteur.fr/listeria).

Results and discussion

A total of 948 tonsil, udder, and feces samples were collected from 316 dairy ruminants (183 cows, 122 sheep, and 11 goats) in two slaughterhouses and one necropsy room over 18 months. Three different Listeria spp. were detected in 35 animals (Table 1 and Table S1). Lm was isolated from 4 (1.3%) and 2 (0.6%) out of 316 animal tonsils, and feces samples, respectively (Table 1 and Table S1).

Liv was isolated from 2 (0.6%), 1 (0.3%), and 13 (4.1%) out of 316 animal tonsils, feces, and udder samples, respectively (Table 1 and Table S1). Lin was isolated from 15 (4.7%) out of 316 feces samples but was not detected in any of the tonsils or udder samples (Table 1 and Table S1). Co-occurrence of different Listeria spp.

(Lm and Lin) was detected in 0.6% (2/316) of individual animal feces (Fig. 1, Fig. S1, Table S1). Lm and Liv were detected more frequently in cattle (19/183 positive cows, prevalence 10.4%) than in small ruminants (2/122 positive sheep, prevalence 1.6%, and 0/11 positive goats, 0%) (ꭕ² test, P<0.05), in agreement with our previous longitudinal study performed in 19 dairy ruminant farms [START_REF] Palacios-Gorba | Ruminant-associated Listeria monocytogenes isolates belong preferentially to dairy-associated hypervirulent clones: a longitudinal study in 19 farms[END_REF].

Regarding, non-pathogenic Listeria spp. the prevalence was significantly higher in small ruminants than in cattle (ꭕ² test, P<0.05), as previously shown [START_REF] Palacios-Gorba | Ruminant-associated Listeria monocytogenes isolates belong preferentially to dairy-associated hypervirulent clones: a longitudinal study in 19 farms[END_REF]. These results show that dairy ruminants can silently carry Lm or Liv in udders or tonsils even without fecal shedding, as we previously reported in wild animals [START_REF] Palacios-Gorba | Listeria spp. Isolated from Tonsils of Wild Deer and Boars: Genomic Characterization[END_REF]. Other studies have shown that large numbers of Lm can be isolated in the tonsils of experimental orally infected sheep despite fecal excretion stopped [START_REF] Zundel | Listeria monocytogenes translocates throughout the digestive tract in asymptomatic sheep[END_REF],

suggesting that tonsils are a location of Lm colonization. Lm is conveyed to milk through contamination of the udder surface or the equipment used for milking [START_REF] Castro | Occurrence, Persistence, and Contamination Routes of Listeria monocytogenes Genotypes on Three Finnish Dairy Cattle Farms: a Longitudinal Study[END_REF].

Lm is infrequently isolated as a cause of dairy ruminant mastitis, which tend to be asymptomatic [START_REF] Watts | Etiological agents of bovine mastitis[END_REF][START_REF] Papić | Source tracking on a dairy farm reveals a high occurrence of subclinical mastitis due to hypervirulent Listeria monocytogenes clonal complexes[END_REF]. The rarity of listerial mastitis cases support the hypothesis that most of the Listeria milk contamination cases may arise from the udder surface or the milking equipment.

Liv positive cows originated more frequently from farms situated in the province of Murcia (14/82, prevalence 17.1%) than in the province of Valencia (2/88, prevalence 2.2%) (P<0.05) (Fig. S2). According to official statistics, the province of Murcia has 77.346 bovines in a surface of 11.313 km², whereas the province of Valencia has 25.664 bovines in a surface of 10.841,82 km² suggesting that Liv prevalence could be influenced by animal density. However, since the number of animals that were slaughtered was not homogeneous across farms (Fig. S2), it could not be determined if the differences in the prevalence of Liv observed between Valencia and Murcia provinces could be attributed to a farm or a province effect. Previous environmental studies have shown that the occurrence of Liv is rare compared with Lm [START_REF] Orsi | Characteristics and distribution of Listeria spp., including Listeria species newly described since 2009[END_REF]. Despite the less frequent isolation of Liv, this species has been found in human feces, domestic ruminant feces and milk, farm environments, dairy processing facilities, and wild boar and deer tonsils [START_REF] Quereda | Pathogenicity and virulence of Listeria monocytogenes: A trip from environmental to medical microbiology[END_REF][START_REF] Hafner | Listeria monocytogenes faecal carriage is common and depends on the gut microbiota[END_REF][START_REF] Alvarez-Ordóñez | Occurrence, Persistence, and Virulence Potential of Listeria ivanovii in Foods and Food Processing Environments in the Republic of Ireland[END_REF][START_REF] Rodriguez | Incidence of Listeria monocytogenes and other Listeria spp. in Ewes' Raw Milk[END_REF]. In line with our results showing the presence of Liv in cow udders, a survey to evaluate the occurrence of Liv in food processing environments and foods in the Republic of Ireland, showed that a higher prevalence was detected for the dairy sector compared to the meat, seafood, and fresh-cut vegetable sectors (1.7 vs 0.2, 0.3, and 0%, respectively) [START_REF] Alvarez-Ordóñez | Occurrence, Persistence, and Virulence Potential of Listeria ivanovii in Foods and Food Processing Environments in the Republic of Ireland[END_REF]. Since a two-step enrichment protocol was used to enhance detection of Listeria spp., and the culled cows did not show mastitis clinical signs in our study, it is plausible that the number of Liv in the analyzed udders was low. The probable low numbers of Liv in the udders together with a restricted geographical distribution of Liv could explain why this species is not commonly detected during surveillance of milk products and only few reports describe its higher prevalence in the milk sector compared to other food sectors [START_REF] Alvarez-Ordóñez | Occurrence, Persistence, and Virulence Potential of Listeria ivanovii in Foods and Food Processing Environments in the Republic of Ireland[END_REF]. Interestingly, Liv isolates obtained from varied geographical and temporal sources in the present study were mainly associated with udders.

Genome-based genotyping results indicated that Lm isolates belonged to CC9 (n=2), CC4 (n=1), CC6 (n=1), CC37 (n=1), and CC388 (n=1) (Table S1, Fig. 1 and Fig. S1). CC9 was previously reported as a CC associated with meat products [START_REF] Maury | Hypervirulent Listeria monocytogenes clones' adaption to mammalian gut accounts for their association with dairy products[END_REF]. CC4 and CC37 were previously reported as CCs associated with dairy products [START_REF] Maury | Hypervirulent Listeria monocytogenes clones' adaption to mammalian gut accounts for their association with dairy products[END_REF]. Importantly, CC4, CC6 and CC388 belong to lineage I which is significantly associated with a clinical origin both in cattle and humans [START_REF] Maury | Uncovering Listeria monocytogenes hypervirulence by harnessing its biodiversity[END_REF].

CC4 and CC6 are more invasive and colonize better the intestinal lumen [START_REF] Maury | Hypervirulent Listeria monocytogenes clones' adaption to mammalian gut accounts for their association with dairy products[END_REF].

Moreover, CC4 and CC6 have been associated with human clinical cases some even without immunosuppressive conditions [START_REF] Maury | Uncovering Listeria monocytogenes hypervirulence by harnessing its biodiversity[END_REF]. This and previous studies suggest that the population of Liv is genetically less diverse than Lm (Fig. S1) [START_REF] Palacios-Gorba | Listeria spp. Isolated from Tonsils of Wild Deer and Boars: Genomic Characterization[END_REF][START_REF] Ramage | Characterisation of Listeria ivanovii isolates from the UK using pulsed-field gel electrophoresis[END_REF]. The significance of these observations and their relation to host adaptation will require further research.

The core virulence genes of LIPI-1 (necessary for the intracellular survival) (the term LIPI stands for Listeria pathogenicity island) and the inlAB locus (necessary for bacterial internalization into eukaryotic cells) were present in all Lm and Liv isolates (Fig. 1). Listeria pathogenic island 3 (LIPI-3), which promotes intestinal colonization by Lm [START_REF] Quereda | Pathogenicity and virulence of Listeria monocytogenes: A trip from environmental to medical microbiology[END_REF], was present in 33.3% (2/6) Lm and 53.3% (8/15) Lin.

Listeria pathogenic island 4 (LIPI-4), which increases Lm neural and placental tropism [START_REF] Maury | Uncovering Listeria monocytogenes hypervirulence by harnessing its biodiversity[END_REF], was present in 33.3% (2/6) Lm and 93.3% (14/15) Lin (Fig. 1). LIPI-2, involved in virulence [START_REF] Quereda | Pathogenicity and virulence of Listeria monocytogenes: A trip from environmental to medical microbiology[END_REF], was present in 100% (16/16) Liv (Fig. 1). Autolysins (aut, aut_IVb), and LPxTG proteins coding genes (inlE, inlG, inlH, inlJ, inlK and vip) related to invasion were present in Lm (Fig. 1, Table S2). Other genes important for teichoic acids biosynthesis, adherence, regulation, and biofilm formation were present in Lm, Liv and Lin (Table S2). Intrinsic antibiotic resistance genes (fosX, Iin (lmo0919), norB, and sul) were present in all Lm, Liv and Lin isolates. tetM and lnuG (resistance to tetracyclines and lincosamides) genes were also present in Lin (5/15, 33.3%; 1/15, 6.7%, respectively). The SSI-1 stress survival islet, involved in Lm growth in high salt concentrations and low pH, was present in 2 out of 3 lineage II Lm isolates. The SSI-2 stress survival islet, which contributes to the alkaline and oxidative stress responses, was present in 14 out of 15 Lin isolates (Fig. 1).

The results of this study contributes to the understanding of the pathogenic potential and ecology of Listeria spp. in dairy ruminants' organs. We have revealed that (i) tonsils and udders constitute a reservoir of Lm or Liv in domestic dairy ruminants; (ii) dairy ruminants can silently carry Lm or Liv in udders or tonsils even without fecal shedding; and (iii) Lm and Liv are detected more frequently in cattle than in small ruminants. Dairy ruminants may contaminate milk or crops and the environment which can pass Lm and Liv to humans via food. The present discoveries could help public health authorities to better control the health risks associated with Lm and Liv.

Fig 1 .

 1 Fig 1. Genomic diversity of Listeria spp. found in this study (N=37). A) Midpoint rooted maximum likelihood phylogeny (GTR+F+I+G4 model, 1000 ultrafast bootstraps) inferred based on the core SNP alignment (327,775 nucleotides) using L. ivanovii subsp. ivanovii CLIP 2021/00402 arbitrary refence. Listeria species are shown in different colors. Information on the isolate's name, origin, date of isolation, CT, and SL are provided. Bold names highlight isolates from the same animal. B) resistance and virulence gene patterns. Color-filled boxes represent the presence of different genetic traits.

  Fig S1. Unrooted core genome phylogenies of Listeria spp. found in this study. A) L. monocytogenes (2,506,493 nucleotides alignment). B) L. ivanovii subsp. ivanovii (2,791,218 nucleotides alignment). C) L. innocua (2,446,398 nucleotides alignment). Tips are colored by source type and sample origin. Dashed circles delimitate clonal complexes (CCs). Isolates obtained from the same animal are marked with a #.

Table 1 .

 1 Liv prevalenceMurcia (14/82 positive cattle, 17.1% Liv prevalence) Number of isolates obtained in this study (N =37).

	Udder
	Feces
	Udder
	Udder
	Udder
	Udder
	Udder
	Udder
	Tonsil
	Udder
	Udder
	Udder
	Udder
	Udder
	Udder
		SL37
		SL9
		SL9
		SL6
		SL387
		SL4
	Feces
	Feces	-
	Feces
	Feces
	Feces
	Feces
	Feces
	Feces
	Feces
	Feces
	Feces
	Feces
	Feces
	Feces
	Feces

Table S1 .

 S1 Listeria isolates characterized in this study.

	BIGSdb id	isolate	aliases	species	isolation date	province	farm	source	source details	no. contigs total length (bp)	N50 contig length (bp)	G+C content (%)	cgMLST loci called (%)	phylogenetic lineage genoserogroup CC (MLST) ST (MLST)	SL (cgMLST)	CT (cgMLST)
	77387	CLIP 2021/00361	LISVAL 368	L. innocua	2019-12-01	Toledo	Farm 16	Animal 8	sheep -feces	26	2893822	579505	37.3	n.a.	n.a.	L	CC493	ST493	n.a.	n.a.
	77388	CLIP 2021/00362	LISVAL 369	L. innocua	2019-12-01	Toledo	Farm 16	Animal 7	sheep -fece	31	2894284	566880	37.3	n.a.	n.a.	L	CC493	ST493	n.a.	n.a.
	77389	CLIP 2021/00363	LISVAL 370	L. innocua	2019-12-01	Toledo	Farm 16	Animal 4	sheep -feces	26	2894491	483337	37.3	n.a.	n.a.	L	CC493	ST493	n.a.	n.a.
	77390	CLIP 2021/00364	LISVAL 371	L. innocua	2019-12-01	Toledo	Farm 16	Animal 6	sheep -feces	31	2894949	270613	37.3	n.a.	n.a.	L	CC493	ST493	n.a.	n.a.
	77391	CLIP 2021/00365	LISVAL 376	L. innocua	2019-12-01	Cuenca	Farm 2	Animal 34	sheep -feces	17	2871184	543767	37.3	n.a.	n.a.	L	CC1480	ST1480	n.a.	n.a.
	77392	CLIP 2021/00366	LISVAL 377	L. innocua	2019-12-01	Cuenca	Farm 3	Animal 39	sheep -feces	20	2784629	479655	37.4	n.a.	n.a.	L	CC599	ST2753	n.a.	n.a.
	77422	CLIP 2021/00378	LISVAL 454	L. innocua	2020-02-13	Murcia	Farm 7	Animal 46	cow -feces	95	2822874	100830	37.4	n.a.	n.a.	L	CC474	ST474	n.a.	n.a.
	77425	CLIP 2021/00381	LISVAL 457	L. innocua	2020-02-13	Cuenca	Farm 8	Animal 67	sheep -feces	67	2871346	144067	37.4	n.a.	n.a.	L	CC530	ST530	n.a.	n.a.
	77426	CLIP 2021/00382	LISVAL 544	L. innocua	2020-02-28	Ciudad Real	Farm 9	Animal 5	goat -feces	46	2847186	140717	37.4	n.a.	n.a.	L	CC530	ST530	n.a.	n.a.
	77427	CLIP 2021/00383	LISVAL 545	L. innocua	2020-02-28	Ciudad Real	Farm 9	Animal 3	goat -feces	45	2846589	188483	37.4	n.a.	n.a.	L	CC530	ST530	n.a.	n.a.
	77428	CLIP 2021/00384	LISVAL 546	L. innocua	2020-02-28	Ciudad Real	Farm 9	Animal 9	goat -feces	48	2846234	219660	37.4	n.a.	n.a.	L	CC530	ST530	n.a.	n.a.
	77432	CLIP 2021/00388	LISVAL 615	L. innocua	2020-06-18	Valencia	Farm 15	Animal 94	cow -feces	28	2971457	333301	37.3	n.a.	n.a.	L	CC2754	ST2754	n.a.	n.a.
	79916	CLIP 2021/01137	LISVAL 772	L. innocua	2021-05-10	Castellón	Farm 13	Animal 119	sheep -feces	38	2815802	138611	37.3	n.a.	n.a.	L	CC1616	ST1616	n.a.	n.a.
	79917	CLIP 2021/01138	LISVAL 773	L. innocua	2021-05-10	Castellón	Farm 13	Animal 120	sheep -feces	46	2864904	210079	37.3	n.a.	n.a.	L	CC1480	ST1480	n.a.	n.a.
	79935	CLIP 2021/01149	LISVAL 784	L. innocua	2021-05-20	Valencia	Farm 14	Animal 183	cow -feces	13	2857362	459131	37.5	n.a.	n.a.	L	CC1482	ST1482	n.a.	n.a.
	77393	CLIP 2021/00367	LISVAL 408	L. monocytogenes	2020-01-16	Murcia	Farm 4	Animal 11	cow -tonsil	35	3000468	477704	37.8	100	II	IIc	CC9	ST9	SL9	CT9007
	77394	CLIP 2021/00368	LISVAL 409	L. monocytogenes	2020-01-16	Murcia	Farm 5	Animal 13	cow -tonsil	84	3003383	192819	37.8	99.6	II	IIc	CC9	ST9	SL9	CT9007
	77395	CLIP 2021/00369	LISVAL 410	L. monocytogenes	2020-01-16	Murcia	Farm 4	Animal 26	cow -tonsil	37	2887721	501910	38	99.9	I	IVb	CC4	ST4	SL4	CT9158
	77421	CLIP 2021/00377	LISVAL 453	L. monocytogenes	2020-02-13	Murcia	Farm 7	Animal 46	cow -feces	67	2991961	150603	37.9	99.9	II	IIa	CC37	ST37	SL37	CT1618
	79286	CLIP 2021/01103	LISVAL 738	L. monocytogenes	2021-04-29	Castellón	Farm 13	Animal 112	sheep -tonsi	43	2926934	292321	37.9	99.9	I	IVb	CC6	ST6	SL6	CT9763
	79918	CLIP 2021/01139	LISVAL 774	L. monocytogenes	2021-05-10	Castellón	Farm 13	Animal 120	sheep -feces	112	2927659	70512	37.9	99.7	I	IVb	CC388	ST388	SL387	CT9761
	77396	CLIP 2021/00370	LISVAL 446	L. ivanovii subsp. ivanovii	2020-02-13	Murcia	Farm 6	Animal 48	cow -tonsil	56	2907510	323421	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77417	CLIP 2021/00373	LISVAL 449	L. ivanovii subsp. ivanovii	2020-02-13	Murcia	Farm 6	Animal 48	cow -udder	54	2926679	348239	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77420	CLIP 2021/00376	LISVAL 452	L. ivanovii subsp. ivanovii	2020-02-13	Murcia	Farm 6	Animal 51	cow -udder	70	2926607	181775	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77423	CLIP 2021/00379	LISVAL 455	L. ivanovii subsp. ivanovii	2020-02-13	Murcia	Farm 6	Animal 57	cow -udder	57	2929111	247159	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77424	CLIP 2021/00380	LISVAL 459	L. ivanovii subsp. ivanovii	2020-02-13	Murcia	Farm 7	Animal 49	cow -udder	55	2913115	352213	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77429	CLIP 2021/00385	LISVAL 605	L. ivanovii subsp. ivanovii	2020-06-02	Valencia	Farm 15	Animal 149	cow -udder	41	2900789	331742	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77433	CLIP 2021/00389	LISVAL 622	L. ivanovii subsp. ivanovii	2020-12-04	Murcia	Farm 1	Animal 65	cow -udder	62	2929665	109996	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77435	CLIP 2021/00391	LISVAL 624	L. ivanovii subsp. ivanovii	2020-12-04	Murcia	Farm 11	Animal 132	cow -udder	48	2901466	164327	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77436	CLIP 2021/00392	LISVAL 625	L. ivanovii subsp. ivanovii	2020-12-04	Murcia	Farm 12	Animal 52	cow -udder	83	2954649	148669	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77437	CLIP 2021/00393	LISVAL 626	L. ivanovii subsp. ivanovii	2020-12-04	Murcia	Farm 10	Animal 145	cow -udder	36	2944244	211155	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77439	CLIP 2021/00395	LISVAL 628	L. ivanovii subsp. ivanovii	2020-12-04	Murcia	Farm 10	Animal 142	cow -udder	41	2901256	224777	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77441	CLIP 2021/00397	LISVAL 630	L. ivanovii subsp. ivanovii	2020-12-04	Murcia	Farm 5	Animal 139	cow -udder	43	2944941	171888	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77443	CLIP 2021/00399	LISVAL 632	L. ivanovii subsp. ivanovii	2020-12-04	Murcia	Farm 10	Animal 148	cow -udder	104	2895532	64194	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77445	CLIP 2021/00401	LISVAL 634	L. ivanovii subsp. ivanovii	2020-12-04	Murcia	Farm 10	Animal 147	cow -tonsil	111	2894126	78390	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	77447	CLIP 2021/00402	LISVAL 635	L. ivanovii subsp. ivanovii	2020-12-04	Murcia	Farm 10	Animal 141	cow -feces	72	2947364	210627	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.
	79911	CLIP 2021/01132	LISVAL 767	L. ivanovii subsp. ivanovii	2021-05-06	Valencia	Farm 15	Animal 162	cow -udder	34	2901547	239098	37	n.a.	n.a.	L	CC883	ST883	n.a.	n.a.

Feces 2020-02-13 2020-02-13 2020-02-13 2020-02-28 2020-02-28 2020-02-28 2020-01-16 2020-01-16 2020-01-16 2020-12-04 2020-12-04 2020-12-04 2020-12-04 2020-12-04 2020-12-04 2019-12-01 2019-12-01 2019-12-01 2019-12-01 2019-12-01 2019-12-01 2020-12-04 2020-12-04 2020-12-04 2021-05-06 2021-04-29 2021-05-10 2021-05-10 2021-05-10 2021-06-18 2021-05-20 2020-06-02 2020-02-13 2020-02-13 2020-02-13 2020-02-13 2020-02-13 Feces

Acknowledgments

This work was supported by Generalitat Valenciana (Project reference AICO/2021/278) (JJQ), the Spanish Ministry of Science and Innovation (Project references PID2019-110764RA-I00/AEI/10.13039/501100011033 (JJQ) and PID2020-119462RA-I00/AEI/10.13039/501100011033 (AGM)) and Universidad CEU Cardenal Herrera Programa INDI 20/40 (JJQ), Institut Pasteur, Inserm, and Santé Publique France (ML). J.J. Quereda is supported by a "Ramón y Cajal" contract of the Spanish Ministry of Science, Innovation, and Universities (RYC-2018-024985-I). C. Palacios-Gorba and Yuval Markovich are supported by a Predoctoral contract from the Universidad Cardenal Herrera-CEU. The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication. We thank the collaboration of Juan M. Corpa who contributed to the sampling of some of the studied sheep and goat in the necropsy room of the Veterinary Faculty.


Guillaume Deltreil
email: guillaume.deltreil@univ-angers.fr

Patrick Tardivel

Piotr Graczyk

Mikael Escobar-Bach

Alexis Descatha
email: alexis.descatha@inserm.fr

EXPOSURE MATRICES WITH JOB HISTORY TO ACCESS WORK EXPOSURE FOR MUSCULOSKELETAL DISORDERS?

Keywords: Occupational, Musculskeletal, pain, lifecourse, mathematical modeling

Introduction. Musculoskeletal disorders related to work might be caused by the cumulative effect of occupational exposures during working life. We aimed to develop a new model which allows to compare the accuracy of duration of work and intensity/frequency associations in application to severe knee pain. Methods. From the CONSTANCES cohort, 66,553 subjects who were working at inclusion and coded were included in the study. The biomechanical job exposure matrix "JEM Constances" was used to assess the intensity/frequency of heavy lifting and kneeling/squatting at work together with work history to characterize the association between occupational exposure and severe knee pain. An innovative model G was developed and evaluated, allowing to compare the accuracy of duration of work and intensity/frequency associations. Results. The mean age was 49 years at inception with 46 percent of women. The G model developed was slightly better than regular models. Among the men subgroup, odds ratios of the highest quartile for the duration and low intensity were not significant for both exposures, whereas intensity/duration were for every duration. Results in women were less interpretable. Conclusion. Though higher duration increased strength of association with severe knee pain, intensity/frequency were important predictors among men. Exposure estimation along working history should have emphasis on such parameters, though other outcomes should be studied and have a focus on women.

Introduction

Musculoskeletal disorders (MSDs) related to repetitive and physically demanding working conditions continue to represent one of the largest occupational disease in industrialized countries [START_REF] Hulshof | The prevalence of occupational exposure to ergonomic risk factors: A systematic review and meta-analysis from the WHO/ILO Joint Estimates of the Work-related Burden of Disease and Injury[END_REF][START_REF] Yassi | Repetitive strain injuries[END_REF]. MSD related to work are caused by non-traumatic injuries, with a possible cumulative effect of occupational exposures during working life, mainly for degenerative disorders like osteoarthritis, and severe knee pain [START_REF] Descatha | Occupational determinants of musculoskeletal disorders. Handbook of disability[END_REF][START_REF] Verbeek | Occupational exposure to knee loading and the risk of osteoarthritis of the knee: a systematic review and a dose-response meta-analysis[END_REF].

The evaluation of biomechanical exposures can be done in different ways and expressed around two main dimensions, intensity/frequency per day and duration over the year [START_REF] Van Der Beek | Assessment of mechanical exposure in ergonomic epidemiology[END_REF]. The exposure assessment can be obtained by estimations based on subjective judgments (self-reports, expert judgments), systematic observations (observations at the workplace, video recording), and direct measurements (at the workplace or in laboratory). However, these techniques are problematic when past exposure evaluations are needed.

Job exposure matrices (JEMs) are commonly used in occupational epidemiology research for the evaluation of past exposures [START_REF] Siemiatycki | Availability of a new job-exposure matrix (CANJEM) for epidemiologic and occupational medicine purposes[END_REF]. Indeed, JEMs allow estimating participants exposures to occupational factors based on job titles, industry sector, and population exposure data. Several biomechanical JEMs have become available recently [START_REF] Vad | Occupational mechanical exposures and reoperation after first-time inguinal hernia repair: a prognosis study in a male cohort[END_REF][START_REF] Rubak | An expert-based job exposure matrix for large scale epidemiologic studies of primary hip and knee osteoarthritis: the Lower Body JEM[END_REF][START_REF] Dalbøge | Upper arm elevation and repetitive shoulder movements: a general population job exposure matrix based on expert ratings and technical measurements[END_REF][START_REF] Evanoff | Using job-title-based physical exposures from O* NET in an epidemiological study of carpal tunnel syndrome[END_REF][START_REF] Dale | General population job exposure matrix applied to a pooled study of prevalent carpal tunnel syndrome[END_REF]. A cumulative exposure index is commonly used to assess cumulative work exposure by multiplying duration and intensity/frequency. However, it is not clear how to consider the combination of intensity/frequency when assessing exposure over the years. Thus, before optimizing models using relevant statistical methods, we first aimed to determine if low level exposure with high duration is equivalent to high level exposure with low duration in the example of severe knee pain and two occupational exposures: carrying heavy loads and kneeling/squatting.

We then aimed develop a new model and compare the accuracy of duration of work and intensity/frequency associations in application to the knee disorders using severe knee pain as an outcome in a large cohort study, by developing an innovative model.

Materials and Methods

Population

The CONSTANCES study is a French general population-based cohort [START_REF] Goldberg | CONSTANCES: a general prospective population-based cohort for occupational and environmental epidemiology: cohort profile[END_REF]. More than 200,000 participants, aged 18-69 years, were recruited between 2012 and 2020 in 23 health screening centers across France. The recruitment was limited to people affiliated to the French National Health Insurance Fund who correspond to active or former salaried workers and their families and excludes agricultural and self-employed workers [START_REF] Goldberg | CONSTANCES: a general prospective population-based cohort for occupational and environmental epidemiology: cohort profile[END_REF]. At enrolment, self-administered questionnaires were sent to participants to collect data including lifestyle, life events, health, and occupations. Variables of interest were collected from the baseline self-administered questionnaire and the medical interview.

For this work, we used French CONSTANCES clean data from 2020. Subjects from this cohort were active at their inclusion with work trajectory coded.

Variables of interest

Participants' sex, age at inception, known inflammatory disease of the joints, regular leisure activity (sports, gardening yes/no), Center for Epidemiologic Studies-Depression Scale (CES-D) into two categories (yes/no), were retrieved from the baseline questionnaire, and body mass index from the medical examination. JEM Constances, which is based on self-reported exposure was used to evaluate occupational exposure [START_REF] Evanoff | The CONSTANCES job exposure matrix based on self-reported exposure to physical risk factors: development and evaluation[END_REF]. In the JEM, occupational exposure is rated from 0 to 4 for intensity/ frequency of heavy lifting ("lifting") and 1 to 4 of kneeling of squatting ("kneeling") based on reported job titles. The JEM Constances was combined with participants reported work trajectory that were coded at baseline retrospectively.

The main outcome was reporting severe knee pain, collected from the self-reported questionnaire at inception: yes if knee pain intensity >5/10 or having knee pain for more than a month per year.

Mathematical modeling

In order to study the influence of duration and intensity/frequency on the onset of the disease, several logistic regression models were built for the two exposures separetely: heavy lifting and kneeling.

The method we propose in this article is based on a generalization of the logistic regression approach. Formally, we define (Y, X i1 , . . . , X ip ) [1,n] an i.i.d sample in {0, 1} × R p of size n ∈ N * where Y is the response variable and corresponds to the illness status of each subject (if the subject is sick Y = 1, if else Y = 0) and X i = (X i1 , . . . , X ip ) are the variables X i1 , . . . , X ip-6 corresponding to a score based on occupation times and levels exposures and X ip-5 , . . . , X ip corresponding to others adjustements variables.

For a given i conditionally on X i := (X i1 , ..., X ip ), Y i follows a Bernoulli distribution such as:

P(Y i = 1|X i ) = π β * (X i ) := exp(< β * , X i >) 1 + exp(< β * , X i >) , (1)
where β * = (β * 0 , β * 1 , . . . , β * p ) ∈ R p is an unknown vector of parameters to be estimated. We estimate the parameter β * , given in [START_REF] Hulshof | The prevalence of occupational exposure to ergonomic risk factors: A systematic review and meta-analysis from the WHO/ILO Joint Estimates of the Work-related Burden of Disease and Injury[END_REF], by Maximal Likelihood Estimation, i.e. by minimizing the normalized opposite of the log likelihood γ n (β) over R p : β = arg min β∈R p γ n (β).

Statistical models

For the occupational health data from the cohort CONSTANCES, several logistic models were possible, depending on the total duration value of the careers and the average exposures of the individuals. We describe them in our cohort via the variables T i = (T i1 , . . . , T ik i ), N i = (N i1 , . . . , N ik i ) and a 1i , . . . , a 6i respectively for occupation times, levels of exposures for each k i jobs held by the subject i as well as six adjustment variables (a 1i = sex, a 2i = age, a 3i = imc, a 4i = leisure, a 5i = arthrite, a 6i = depression). We must take into account that the number of jobs k i can be very different. The models A,B and C are defined through the three following transformations:

X A i = m ∑ j=1 N ij T ij X B i = m ∑ j=1 N B ij T ij , with N B ij = 0, if N ij ≤ 1 N ij , if else. X C i = m ∑ j=1 N C ij T ij , with N C ij = 0, if N ij ≤ 3 1, if else. Thus, Y i follows a Bernoulli distribution such that P(Y i = 1|X i ) = logit(β * 0 + β * 1 X • i + ∑ 6 j=1 β * j+1 a ji ) where X • i is to be replaced by X A i , X B i or X C
i depending on the variable considered. We specify that the model B differs, from the model A, because in the computation of its transformation the exposure levels between 0 and 1 are confounded. The model C differs, from the model A, because in the computation of its transformation the lowest exposure levels are nullified. In order to define the model G we need to introduce two transformations ϕ 1 and ϕ 2 :

• ϕ 1 (T i ) = ∑ k i
j=1 T ij corresponding to the total duration of the career range.

• ϕ 2 (T i , N i ) = ∑ k i j=1 N ij g(T ij ) ∑ k i j=1 g(T ij )
, where g(x) = 1exp (-x) which is the time-weighted average of the exposure level.

The construction of the design matrix (X ij ) (1≤i≤n,1≤j≤p) of the model G is based on a semi-coding of the variables {ϕ 1 (T i )} 1≤i≤n and {ϕ 2 (T i , N i )} 1≤i≤n . More precisely, the membership in a group is determined by the belonging of the values ϕ 1 (T i ) and ϕ 2 (T i , N i ) to different given intervals. For this study, we consider the classes of intervals

G 1 := {[q j , q j+1 [, j = 1, . . . , q} and G 2,c := { [j -1, j[ } c≤j≤4 ∪ {4}, c = 1, 2
where q j ∈ N * defines fixed empirical quantiles of the sample {ϕ 1 (T i )} 1≤i≤n . When c = 1, the exposure considered is carrying a heavy load, and when c = 2 when the exposure considered is kneeling. The set of groups is defined as the class G := {G j , j = 1, . . . , p} = G 1 × G 2 where p = (6c)q. An individual is then associated with the group

j ∈ [[1, p]] with X ij = 1 if and only if (ϕ 1 (T i ), ϕ 2 (T i , N i )) ∈ G j .

Selection of the design matrix for model G

We propose to choose the design matrix X for the model G by a recent model selection procedure, introduced and described in the article "Model selection in logistic regression" by Kwemou et al. [START_REF] Kwemou | Model selection in logistic regression[END_REF]. The mathematical guarantees for this model selection method are based on oracle inequalities from Birgé and Massart [START_REF] Birgé | Gaussian model selection[END_REF]. This model selection is performed using penalized maximum likelihood estimators which will allow us to choose the best design matrix.

Let F be a family of design matrices X (m) , m ∈ {1, . . . , M}. For each m ∈ {1, . . . , M}, we define βm as the estimator obtained by minimizing γ n over R p m where p m is the number of columns of X (m) , namely: βm = arg min β∈R pm γ n (β) We use a data driven strategy that selects the best matrix among the family F . For this purpose we use a penalized maximum likelihood criterion for choosing the index m associated to the appropriate design. Here, we consider the Akaike information criterion where p m corresponds to the number of parameters in the model (1), i.e. p m is the number of columns of the matrix X (m) and m = arg min m∈{1,...,M}

γ n ( β m ) + p m .
Hence, minimizing this criterion allows us to find the best design matrix X m .

The "lifting" and "kneeling" quartiles of high exposure is considered in terms of duration and intensity/frequency:

• into 5 categories for heavy lifting G 2,1 := { [j -1, j[ } 1≤j≤4 ∪ {4}, • into 4 categories for kneeling G 2,2 := { [j -1, j[ } 2≤j≤4 ∪ {4}.

α-divergence and Statistical tests

The logistic regression formula given in (1) depends on the model A, B, C and G. We wish to compare these statistical models, via a differentiation criterion, in order to identify the most appropriate model for this formula. Let r be the unknown probability that a randomly selected subject from the population is sick. It is natural to estimate the unknown parameter r by the proportion of patients observed in our dataset: r1 = 1 n ∑ n i=1 Y i . One can also estimate r from the logistic regression formula via the estimator r2 = 1 n ∑ n i=1 π β * (X i ). Note that unlike r1 , the estimator r2 depends on the model. For an appropriate model it is natural to choose the model for which the estimators r1 and r2 are close.

We here chose to rely on the work of A.Basu et al. [START_REF] Basu | Robust and efficient estimation by minimising a density power divergence[END_REF] who developped a robust estimator for the density function. The criterion used in this article is based on the αdivergence between two densities f and g (relative to a measure µ) defined for an α > 0 as follows: Here, we need to use the α-divergence between two Bernoulli laws with parameters r 1 ∈]0.1[ and r 2 ∈]0.1[ which we define by

d α (g, f ) = f 1+α (z) -1 + 1 α g(z) f α (z) + 1 α g 1+α (
d α (r 1 , r 2 ) =        r 1+α 2 + (1 -r 2 ) 1+α -1 + 1 α ((1 -r 2 ) α (1 -r 1 ) + r α 2 r 1 ) + 1 α ((1 -r 1 ) 1+α + r 1+α 1 ) α > 0 r 1 ln r 1 r 2 + (1 -r 2 ) ln 1-r 1 1-r 2 α = 0
The parameter α here defines a trade-off between the estimation efficiency and the variability robustness. Note that when α is close to 0 we retrieve the Kullback-Leibler divergence and when α = 1 we have

d 1 (r 1 , r 2 ) = 2(r 1 -r 2 ) 2 .
Next, we compared the p-value of the Wald's test (with Bonferroni correction) and odds ratios (OR) of the highest quartile for the duration with low intensity/frequency and the highest quartile for the intensity/frequency with low duration. Stratification on sex was performed as sensitivity analysis. For both exposures, logistic models were built adjusted on relevant variables. We compared the p-value of the Wald's test (with Bonferroni adjustments) and odds ratios of different quartile of duration and intensity/frequency and with low duration and lowest intensity as reference.

Analysis plan

After a brief description of the sample and the available adjustement variables, we assessed the performance score based on the α-divergence allowing the comparison of different models with α varying with each relevant occupational exposure. For a small value of parameter α, the decision of the model choice can be questioned but when α increases the choice of model becomes easier and the model G is selected.

Then, we were able to compare Odds ratios with the Wald's test. Sex stratification has been performed as primary analysis (whereas 45 year of age stratification has been considered in secondary analysis).

Analyses were performed on Python 3.8.8 (statsmodels,statistics,pandas,numpy).

Results

The sample included 66,553 subjects (Table 1), with 30,765 women (46.2%) and with a median of 49 years old (21 years of employment). The subjects reported regular primary leisure activities in 43.1% of cases, with half of sample who were overweight/ obese, and an important part (21.2%). who had a positive CES-D, suggesting a depression. In order to study intensity and duration of exposure, we first selected the best model between A to G. We compared this familly of models by increasing α. For any value of parameter alpha, the model G is selected (Tables 2 and3).

Figures, Tables and Schemes

For heavy lifting as well as kneeling, intensity/ frequency was the most important predictor. Duration increased risks only for men with a dose response relationship. For women, only intensity and frequency in low duration for both exposures seemed associated with the knee pain (Tables 4 and5). Although the statistical power was lower, the pattern of association was similar for participants aged less and more than 45 years (supplemental Tables 6 and7).

Discussion

This is the first study that used a developed mathematical model to compare the effect of duration and intensity during working life, on the association with severe knee as a proxy of degenerative musculoskeletal disorders. The new model G was found to be better than the usual ones, though the difference was minor. For men, we found that the OR of the highest quartile for the duration and low intensity is not significant for both exposures, whereas intensity/duration is significant for every duration, with a dose response relationship. Results for women were limited.

As expected, there was as an important effect of the intensity of heavy lifting and kneeling on severe knee pain. Both exposures are known to be associated with knee disorders [START_REF] Klussmann | Individual and occupational risk factors for knee osteoarthritis: results of a case-control study in Germany[END_REF][START_REF] Descatha | Long-term effects of biomechanical exposure on severe knee pain in the Gazel cohort[END_REF][START_REF] Fransen | Hip and knee pain: Role of occupational factors[END_REF]. The dose response relationship has been described previously. Jensen calculated an equivalent of our model A using an individual exposure from the number of knee-straining activities and the number of years in the trade within a collective of floor layers, carpenters and compositors. The ORs for knee complaints and radiographically determined knee osteoarthritis were 3.0 (95% CI, 0.5 to 17.2) in the low-exposure group, 4.2 (95% CI, 0.6 to 27.6) in the medium-exposure group, and 4.9 (95% CI, 1.1 to 21.9) in the high-exposure group compared with the zero-exposure group [START_REF] Jensen | Knee-straining work activities, self-reported knee disorders and radiographically determined knee osteoarthritis[END_REF]. There is an important difference in the strength of the associations compared to our work but it should be explained by the large population design with JEM exposure methods. Indeed, high exposure is considered using the proxy of job title with a large variation inside job categories.

Also expected, in a previous review on occupational exposure and knee osteoarthritis, [START_REF] Jensen | Knee osteoarthritis: influence of work involving heavy lifting, kneeling, climbing stairs or ladders, or kneeling/squatting combined with heavy lifting[END_REF] lifting and carrying of loads was significantly associated with severe knee pain. Knee osteoarthritis was also found associated with lifting and carrying of loads with a doseresponse relationship : OR of 2.0 (95% CI, 1.1 to 3.6) in the exposure group 630 to <5,120 kg-hours over life, up to an OR of 2.6 (95% CI, 1.1 to 6.1) in the highest exposure group (>37,000 kg-hours over life) in men [START_REF] Jensen | Knee-straining work activities, self-reported knee disorders and radiographically determined knee osteoarthritis[END_REF].

The lack of clear association for women was also found by D'Souza et al. who reported on an analysis of the US national survey, where they describe relationships between work activities and symptomatic knee osteoarthritis [START_REF] Souza | Analysis of the Third National Health and Nutrition Examination Survey (NHANES III) using expert ratings of job categories[END_REF]. A significant exposure-response relationship was only found between symptomatic knee osteoarthritis and kneeling in men but not women. Different explanations might be suggested: since our model included adjusting factors like BMI and depression, there might be a more complex causal pathways than in men such as considered in back pain [START_REF] Leclerc | Level of education and back pain in France: the role of demographic, lifestyle and physical work factors[END_REF]. JEM Constances is not gender stratified and applying a specific JEM for sex could be a lead for another study. Furthermore, selection effect similar to healthy workers effect is also possible. The main strength of our study was the possibility to use JEM with working life course on a large cohort study. Limitations might also be raised by using a large but not representative population of French workers, and specific jobs in agriculture (not included by design) and mining (almost disappeared in France) should be considered since they are known factors related to lower limb MSD [START_REF] Fransen | Hip and knee pain: Role of occupational factors[END_REF]. Second, the outcome was focused on severe knee pain. Even, it is self-declared and might correspond to heterogenous disorders, work-related or not. It was used as a good example of a proxy of degenerative musculoskeletal disorders, and the use of pain intensity and severity is recommended [START_REF] Hagberg | Prevention of musculoskeletal disorders in workers: classification and health surveillance-statements of the Scientific Committee on Musculoskeletal Disorders of the International Commission on Occupational Health[END_REF]. We have previously shown that working in a kneeling or squatting position was significantly associated with severe knee pain [START_REF] Descatha | Long-term effects of biomechanical exposure on severe knee pain in the Gazel cohort[END_REF]. More recently, similar trends of associations between severe knee pain and knee arthroplasty groups were showed in the same cohort [START_REF] Valter | Influence of severe knee pain, meniscus surgery and knee arthroplasty on physical ability: an observational study of 114 949 adults in the CONSTANCES cohort[END_REF]. This result is also found elsewhere, with non-managerial jobs associated with higher prevalence of knee osteoarthritis and knee symptoms [START_REF] Lee | Effects of education, income, and occupation on prevalence and symptoms of knee osteoarthritis[END_REF]. Third, as we already mentioned, the use of JEMs might also be questioned since it is a global average evaluation that does not consider the differences inside job [START_REF] Descatha | Job-Exposure Matrix: A Useful Tool for Incorporating Workplace Exposure Data Into Population Health Research and Practice[END_REF]. However, assessing exposure during long periods of time and for a big number of subjects is challenging and JEMs are appropriate tool to consider. Furthermore, assessement of carrying heavy loads exposure using JEMs was found to be valid compared to a self-administrated questionnaire [START_REF] Ngabirano | Comparison between a job-exposure matrix (JEM) score and self-reported exposures for carrying heavy loads over the working lifetime in the Constances cohort[END_REF].

Conclusions

This innovative approach using mathematic modeling of working history and a JEM, shows that duration in years has a smaller impact than frequency/intensity and should be considered at least among men. Our new model G seems to be an interesting approach though the improvement is light. Further study, should be done on other outcomes, and have a focus on women.

  z) dµ(z) Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 30 September 2022 doi:10.20944/preprints202209.0489.v1

Table 1 .

 1 Description of the sample.

	Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 30 September 2022	doi:10.20944/preprints202209.0489.v1
	Variables		N(%)	Mean (SD)
	Sex	Men	35788 (53.7)
		Women	30765 (46.2)
	Age (years)			48.5 (13.1)
	Body Mass Index (kg/m²)			25.07 (24.0)
	Leisure activites	Yes	28693 (43.1)
		No	37860 (56.9)
	Inflammatory Ostearthritis	Yes	950 (1.4)
		No	65603 (98.6)
	Depression	Yes	14095 (21.2)
		No	52458 (78.8)

Table 2 .

 2 Comparison of performance score for severe knee pain with lifting.

	α	A	B	C	G
	0	0.305	0.329	0.326	0.192
	0.25	0.162	0.228	0.218	0.019
	0.5	0.155	0.221	0.211	0.017
	0.75	0.146	0.210	0.200	0.015
	1	0.136	0.198	0.189	0.013

Table 3 .

 3 Comparison of performance score for severe knee pain with kneeling.

	α	A	B	C	G
	0	0.282	0.304	0.297	0.194
	0.25	0.101	0.160	0.145	0.021
	0.5	0.096	0.153	0.138	0.018
	0.75	0.092	0.143	0.13	0.016
	1	0.088	0.134	0.121	0.014

Table 4 .

 4 Results of adjusted logistic regression of severe knee pain with lifting, for men and women separately.

	Variable		Men		Women
	Duration	Intensity/Frequency OR	IC 95%	p-value	OR	IC 95%	p-value
	low	[1,2[	1.83	[1.53 , 2.19]	< 10 -4	1.46	[1.29 , 1.66]	< 10 -4
	low	[2,3[	2.08	[1.71 , 2.52]	< 10 -4	1.47	[1.20 , 1.81]	< 10 -4
	low	[3,4[	1.96	[1.42 , 2.70]	< 10 -4	1.38	[ 0.91 , 2.09 ]	0.13
	low	{4}	2.49	[1.86 , 3.33]	< 10 -4	2.81	[1.60 , 4.92]	< 10 -4
	medium	[0,1[	1.07	[0.91 , 1.27]	0.41	0.86	[0.77 , 0.96]	0.01
	medium	[1,2[	1.63	[1.35 , 1.97]	< 10 -4	1.11	[0.97 , 1.28]	0.13
	medium	[2,3[	2.39	[1.97 , 2.90]	< 10 -4	1.28	[1.02 , 1.61]	0.04
	medium	[3,4[	2.45	[1.78 , 3.39]	< 10 -4	1.08	[ 0.62 , 1.86 ]	0.79
	medium	{4}	2.71	[1.84 , 3.98]	< 10 -4	1.09	[0.41 , 2.93]	0.86
	high	[0,1[	1.09	[0.92 , 1.29]	0.33	0.87	[ 0.77 , 0.98]	0.03
	high	[1,2[	1.77	[1.49 , 2.11]	< 10 -4	1.07	[0.93 , 1.23]	0.36
	high	[2,3[	2.12	[1.77 , 2.53]	< 10 -4	1.09	[0.86 , 1.37]	0.49
	high	[3,4[	2.41	[1.87 , 3.11]	< 10 -4	1.37	[0.99 , 1.91]	0.06
	high	{4}	3.56	[2.68 , 4.72]	< 10 -4	1.67	[0.68 , 4.12]	0.27
			Adjusted on body mass index, leisure activity, inflammatory osteoarthritis, depression, and age/sex when not
			stratified ; Reference duration low, Intensity/Frequency [0,1[.

Preprints (www.

preprints.org) | NOT PEER-REVIEWED | Posted: 30 September 2022 doi:10.20944/preprints202209.0489.v1

Table 5 .

 5 Results of adjusted logistic regression of severe knee pain with kneeling, men and women separately.

	Variable		Men			Women
	Duration	Intensity/Frequency OR	IC 95%	p-value	OR	IC 95%	p-value
	low	[2,3[	1.97	[1.67 , 2.31]	< 10 -4	1.34	[1.16 , 1.53]	< 10 -4
	low	[3,4[	1.75	[1.40 , 2.20]	< 10 -4	1.44	[1.22 , 1.71]	<10 -4
	low	{4}	1.98	[1.51 , 2.60]	< 10 -4	1.35	[1.01 , 1.82]	0.04
	medium	[1,2[	1.06	[0.91 , 1.23]	0.48	0.83	[0.74 , 0.93]	<10 -4
	medium	[2,3[	1.86	[1.57 , 2.20]	< 10 -4	1.02	[0.88 , 1.17]	0.83
	medium	[3,4[	2.19	[1.72 , 2.79]	< 10 -4	1.23	[1.03 , 1.48]	0.02
	medium	{4}	2.66	[1.91 , 3.71]	< 10 -4	0.91	[0.59 , 1.39]	0.66
	high	[1,2[	1.10	[0.94 , 1.28]	0.26	0.85	[0.76 , 0.96]	0.01
	high	[2,3[	1.80	[1.54 , 2.11]	< 10 -4	1.03	[0.89 , 1.19]	0.67
	high	[3,4[	2.16	[1.78 , 2.62]	< 10 -4	0.88	[0.74 , 1.04]	0.13
	high	{4}	2.55	[1.99 , 3.28]	< 10 -4	1.38	[1.03 , 1.84]	0.03

Adjusted on body mass index, leisure activity, inflammatory osteoarthritis, depression, and age/sex when not stratified ; Reference duration low, Intensity/Frequency

[1,2[.

Table 6 .

 6 Results of adjusted logistic regression of severe knee pain with lifting, for <45 years old and 45 or more participants.Adjusted on body mass index, leisure activity, inflammatory osteoarthritis, depression, and age/sex when not stratified ; Reference duration low, Intensity/Frequency [0,1[.

	Variable		<45 years			≥45 years
	Duration	Intensity/Frequency OR	IC 95%	p-value	OR	IC 95%	p-value
	low	[1,2[	1.59	[1.35 , 1.87]	< 10 -4	1.51	[1.26 , 1.81]	< 10 -4
	low	[2,3[	1.71	[1.37 , 2.14]	< 10 -4	1.62	[1.31 , 2.00]	< 10 -4
	low	[3,4[	1.45	[0.95 , 2.19]	0.08	1.42	[0.93 , 2.16]	0.107
	low	{4}	2.63	[1.83 , 3.79]	< 10 -4	1.31	[0.86 , 2.02]	0.212
	medium	[0,1[	0.87	[0.75 , 1.02]	0.09	0.92	[0.79 , 1.06]	0.234
	medium	[1,2[	1.33	[1.10 , 1.59]	0.002	1.24	[1.05 , 1.46]	0.01
	medium	[2,3[	1.64	[1.28 , 2.10]	< 10 -4	1.57	[1.28 , 1.93]	< 10 -4
	medium	[3,4[	1.79	[1.19, 2.70]	0.005	1.95	[1.34 , 2.84]	0.001
	medium	{4}	3.30	[2.04, 5.34]	< 10 -4	2.25	[1.42, 3.58]	0.001
	high	[0,1[	0.93	[0.77, 1.12]	0.442	0.88	[0.78, 1.01]	0.073
	high	[1,2[	1.31	[1.07, 1.62]	0.01	1.21	[1.05, 1.39]	0.007
	high	[2,3[	1.93	[1.51, 2.47]	< 10 -4	1.46	[1.25, 1.71]	< 10 -4
	high	[3,4[	1.36	[0.85, 2.19]	0.204	1.69	[1.36, 2.09]	< 10 -4
	high	{4}	1.63	[0.89, 2.98]	0.114	2.58	[1.97, 3.37]	< 10 -4

Table 7 .

 7 Results of adjusted logistic regression of severe knee pain with kneeling, for <45 years old and 45 or more participants.Adjusted on body mass index, leisure activity, inflammatory osteoarthritis, depression, and age/sex when not stratified ; Reference duration low, Intensity/Frequency[1,2[.

	Variable		<45 years			≥45 years
	Duration	Intensity/Frequency OR	IC 95%	p-value	OR	IC 95%	p-value
	low	[2,3[	1.50	[1.28 , 1.78]	< 10 -4	1.49	[1.24 , 1.79]	< 10 -4
	low	[3,4[	1.62	[1.31 , 2.01]	< 10 -4	1.45	[1.16 , 1.80]	0.001
	low	{4}	1.80	[1.33 , 2.44]	< 10 -4	1.40	[1.04, 1.88]	0.025
	medium	[1,2[	0.86	[0.74 , 1.00]	0.06	0.91	[0.79 , 1.05]	0.19
	medium	[2,3[	1.45	[1.22 , 1.73]	< 10 -4	1.20	[1.02 , 1.41]	0.028
	medium	[3,4[	1.34	[1.04 , 1.73]	0.026	1.47	[1.20 , 2.47]	< 10 -4
	medium	{4}	1.24	[0.78, 1.97]	0.361	1.75	[1.24, 2.47]	0.001
	high	[1,2[	0.91	[0.76, 1.09]	0.300	0.88	[0.78, 1.00]	0.050
	high	[2,3[	1.38	[1.14, 1.69]	0.001	1.23	[1.08, 1.41]	0.002
	high	[3,4[	1.59	[1.22, 2.09]	0.001	1.22	[1.05, 1.42]	0.011
	high	{4}	1.42	[0.89, 2.26]	0.139	1.73	[1.41, 2.12]	< 10 -4

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 30 September 2022

Acknowledgments:

We thank Sabrina Pitet for her help in improving manuscript tables and reference list. The authors also thank the team of the "Population-based Epidemiologic Cohorts Unit" (Cohortes en population) that designed and manages the CONSTANCES Cohort Study, as well as the volunteers.

Data Availability Statement: Personal health data underlying the findings of our study are not publicly available due to legal reasons related to data privacy protection. However, the data are available upon request to all interested researchers after authorization of the French "Commission nationale de l'informatique et des libertés". The CONSTANCES email address is contact@constances.fr.

Informed Consent Statement: Informed consent was obtained from all subjects involved in the study

Conflicts of Interest:

The authors declare no conflict of interest (paid by their affiliations. Authors are paid by their institution, and A.D. is also paid as editor of the Archives des Maladies professionnelles et de l'Environnement (Elsevier).The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results".

Abbreviations

The following abbreviations are used in this manuscript: OR odds ratio


Florian Miralles
email: florian.miralles@umontpellier.fr

Bastien Sauvage
email: bastien.sauvage@inria.fr

Stephen Wornom
email: stephen.wornom@inria.fr

Bruno Koobus
email: bruno.koobus@umontpellier.fr

Alain Dervieux
email: alain.dervieux@inria.fr

Application of hybrid RANS/VMS modeling to rotating machines

Keywords: hybrid turbulence models, variational multiscale model for LES, chimera method, cylinder flow, rotating machine, unstructured grids

The proposed communication deals with hybrid RANS-LES modeling. The target application is the study of flows around rotating machines like helicopters and drones. In fine, the simulations should provide accurate estimates concerning the noise emission. Each of these flows can involve mean and high Reynolds turbulent regions with detached eddies and with thin laminar and turbulent boundary layers. A hybrid model, like DDES, is then mandatory, with possibly an improved resolution of LES regions, which are mainly turbulent wakes. It is then interesting to apply there a more sophisticated LES model than the LES part of DDES. In our study, we use there the Dynamic variational multiscale model (DVMS). In the other regions, a DDES or simply a RANS modeling is applied. In both cases a twoequation closure is chosen. After a discussion of the modeling ingredients, we shall present a comparison of the RANS, LES, and hybrid models for two series of flows. Although computed by many researchers, flows around cylinders remain difficult to predict. The comparison will continue with a flow around a cross shaped mixing device rotating inside a cylinder.

INTRODUCTION

In this paper, hybrid models are evaluated for the simulation of massively separated flows around fixed and moving geometries, with the objective of applying them to aeroacoustic problems involving complex industrial flows at high Reynolds numbers. For this purpose, the turbulence models need to behave preferentially like large eddy simulation models (LES) in order to introduce little dissipation and thus better capture the small scales and the fluctuations of the unsteady flows considered. However, these flows involve very thin boundary layers, which current computers and softwares cannot compute in LES mode, only in RANS mode. To accommodate both needs, LES and RANS, the hybrid turbulence approach, combining LES and RANS in a somewhat zonal manner, is considered by many research teams.

In this work, beside the classical Delayed Detached Eddy Simulation model (DDES, [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF]), we study an approach which combines zonally this last model and the dynamic variationnal multiscale model (DVMS, [START_REF] Koobus | A variational multiscale method for the large eddy simulation of compressible turbulent flows on unstructured meshes-application to vortex shedding[END_REF]). In this hybrid DDES/DVMS approach [START_REF] Moussaed | Simulation of the flow past a circular cylinder in the supercritical regime by blending RANS and variational-multiscale LES models[END_REF], the DVMS model is preferentially activated in the wake regions where this latter model introduces less dissipation than the LES component of DDES. A hybrid RANS/DVMS model [START_REF] Itam | Combining a DDES model with a dynamic variational multiscale formulation[END_REF], where the RANS component is Goldberg's k -ε model [START_REF] Goldberg | A wall-distance-free k -ε model with Enhanced Near-Wall Treatment[END_REF], is also applied in this study. A smooth blending function, which is based on the value of a blending parameter, is employed in these hybrid strategies, for switching from RANS to DVMS or DDES to DVMS. In [START_REF] Itam | Combining a DDES model with a dynamic variational multiscale formulation[END_REF], these models are applied to flows around a cylinder at moderate Reynolds number, starting from Reynolds number 20000. In many industrial cases, and particularly cases with rotating geometries like propellers, the Reynolds number is much higher.

In the present paper, these hybrid models are applied to the flow around a circular cylinder in the supercritical regime, namely at Reynolds numbers 1 million and 2 million. At Reynolds number 1 million, the supercritical regime shows turbulent flow separation on both the upper and lower surface of the cylinder with a laminar-turbulent transition in the boundary layers located between the stagnation point and the separation one. At Reynolds number 2 million, the high transition regime shows a fully turbulent boundary layer on one side of the cylinder, and a partially laminar, partially turbulent boundary layer on the other side. This benchmark, which contains many characteristics encountered in industrial flows, is challenging due to the complex physics of the flow and the considered high Reynolds numbers, and remains a difficult calculation to perform. The second application concerns the flow at Reynolds number 1.8 million around a cross shaped mixing device rotating inside a cylinder.

This work is part of a cooperation program with Keldysh Institute of Applied Mathematics of Moscow focussing on the "Efficient simulation of noise of rotating machines", see e.g. [START_REF] Garanzha | Moving deforming mesh generation based on quasi-isometric functional[END_REF][START_REF] Bakhvalov | On using artificial viscosity in edge-based schemes on unstructured grids[END_REF][START_REF] Duben | EBR schemes with curvilinear reconstructions of variables in the near-wall region[END_REF].

The remainder of this paper is organized as follows. Section 2 presents the hybrid turbulence models used in this work. Section 3 describes the numerical scheme. In Section 4, applications are presented. The obtained results are analyzed and compared with those obtained in other numerical studies and with experimental data. Finally, concluding remarks are drawn in Section 4.

TURBULENCE MODELING 2.1. RANS

First, we want to specify that RANS stands for unsteady RANS throughout the document.

In this work, and as far as the closure of the RANS equations is concerned, we use the low Reyn-olds k -ε model proposed in Goldberg et al. [START_REF] Goldberg | A wall-distance-free k -ε model with Enhanced Near-Wall Treatment[END_REF]. This model was designed to improve the predictions of the standard k-ε one for adverse pressure gradient flows, including separated flows. For the sake of brevity, the equations of this model are not recalled in this document.

LES-like model: DVMS

In this section, we present the DVMS model which is preferred to the classical LES approach in our hybrid strategies because of some specific properties, as will be specified hereafter, that allow this model to have a better behavior in some regions of the flow, such as shear layers and wakes.

The Variational Multiscale (VMS) model for the large eddy simulation of turbulent flows has been introduced in [START_REF] Hughes | Large eddy simulation and the variational multiscale method[END_REF] in combination with spectral methods. In [START_REF] Koobus | A variational multiscale method for the large eddy simulation of compressible turbulent flows on unstructured meshes-application to vortex shedding[END_REF], an extension to unstructured finite volumes is defined. This method is adapted in the present work. Let us explain this VMS approach in a simplified context. Assume the mesh is made of two embedded meshes. On the fine mesh we have a P 1 -continuous finite-element approximation space V h with the usual basis functions Φ i vanishing on all vertices but vertex i. Let V 2h represents its embedded coarse subspace. Let V h be the complementary space: V h = V 2h ⊕ V h . The space of small scales V h is spanned by only the fine basis functions Φ i related to vertices which are not vertices of V 2h . Let us denote the compressible Navier-Stokes equations by:

∂W ∂t + ∇ • F(W) = 0 where W = (ρ, ρu, E).
The VMS discretization writes for W h = W i Φ i :

∂W h ∂t , Φ i + ∇ • F(W h ), Φ i = -τ LES (W h ), Φ i (1)
For a test function related to a vertex of V 2h , the RHS vanishes, which limits the action of the LES term to small scales. In practice, embedding two unstructured meshes V h and V 2h is a constraint that we want to avoid. The coarse level is then built from the agglomeration of vertices/cells as sketched in Fig. 1.

It remains to define the model term τ LES (W h ). This term represents the subgrid-scale (SGS) stress term, acting only on small scales W h , and computed from the small scale component of the flow field by applying either a Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF] or a WALE SGS model [START_REF] Nicoud | Subgridscale stress modelling based on the square of the velocity gradient tensor[END_REF], the constants of these models being evaluated by the Germano-Lilly dynamic procedure [START_REF] Germano | A dynamic subgridscale eddy viscosity model[END_REF][START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF].

The resulting model is denoted DVMS in this paper. It has been checked [START_REF] Moussaed | Impact of dynamic subgrid-scale modeling in variational multiscale large-eddy simulation of bluff body flows[END_REF] that combining VMS and dynamic procedure effectively brings improved predictions.

A key property of the VMS formulation is that the modeling of the dissipative effects of the unresolved structures is only applied on the small resolved scales, as sketched in Fig. 2. This property is not satisfied by LES models which damp also the other resolved scales. Important consequences are that a VMS model introduces less dissipation than its LES counterpart (based on the same SGS model) and that the backscatter transfer of energy from smallest scales to large scales is not damped by the model. These VMS models then generally allow a better behavior near walls, in shear layers and in the presence of large coherent structures. Moreover, in this work, the dynamic procedure, which provides a tuning of the SGS dissipation in space and time, is combined with the VMS approach, which limits its effects to the smallest resolved scales, so that the resulting DVMS model enjoys synergistic effects.

DDES model

The DDES approach used in this work is based on the Goldberg's k-ε model [START_REF] Goldberg | A wall-distance-free k -ε model with Enhanced Near-Wall Treatment[END_REF]. In other words, this RANS model is introduced in the DDES formulation [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF] l DDES where

l DDES = k 3/2 ε -f d max(0, k 3/2 ε -C DDES ∆) with f d = 1 -tanh((8r d ) 3 )
and

r d = ν t + ν max( √ u i, j u i, j , 10 -10 )K 2 d 2 w .
K denotes the von Karman constant (K = 0.41), d w the wall-normal distance, u i, j the x j -derivative of the ith-component of the velocity u, and the model constant C DDES is set to the standard value 0.65.

Hybrid RANS/DVMS

The central idea of the hybrid RANS/DVMS model [START_REF] Moussaed | Simulation of the flow past a circular cylinder in the supercritical regime by blending RANS and variational-multiscale LES models[END_REF] applied in this study is to combine the mean flow field obtained by the RANS component with the application of the DVMS model wherever the grid resolution is adequate. In this hybrid model, the k -ε model of Goldberg (subsection 2.1) is used as the RANS component. First, let us write the semidiscretization of the RANS equations :

∂ W h ∂t , Φ i + ∇ • F( W h ), Φ i = -τ RANS ( W h ), Φ i .
A natural hybridation writes:

∂W h ∂t , Φ i + ∇ • F(W h ), Φ i = -θ τ RANS ( W h ), Φ i -(1 -θ) τ LES (W h ), Φ i
where W h denotes now the hybrid variables and θ is the blending function which varies between 0 and 1 and is defined by:

θ = 1 -f d (1 -tanh(ξ 2 )) with ξ = ∆ l RANS or ξ = µ SGS µ RANS ,
the shielding function f d being defined as in Subsection 2.3.

The blending function θ allows for a progressive switch from RANS to LES (DVMS in our case) where the grid resolution becomes fine enough to resolve a significant part of the local turbulence scales or fluctuations, i.e. computational regions suitable for LES-like simulations. Additionally, this blending function, thanks to the shielding function f d , prevents the activation of the LES mode in the boundary layer. By way of example, the isocontours of the blending function θ for the flow past a circular cylinder are shown in Fig. 3. We can see in particular that in the boundary layer the RANS model is activated, while in the wake the LES approach is recovered. for θ = 1 (in red) the RANS model is activated, wherever 0 < θ < 1 additional resolved fluctuations are computed, and in the limit θ → 0 (in blue) the LES approach (DVMS in this work) is recovered.

Hybrid DDES/DVMS

The key idea of the proposed hybrid DDES/DVMS model [START_REF] Itam | Combining a DDES model with a dynamic variational multiscale formulation[END_REF] is to use the DVMS approach instead of the LES component of DDES in locations where this component is expected to be activated, especially in wake regions where the DVMS approach allows more accurate prediction.

Assuming that the semi-discretization of the DDES equations writes:

∂W h ∂t , Φ i + ∇ • F(W h ), Φ i = -τ DDES (W h ), Φ i ,
the hybrid equations are then defined by:

∂W h ∂t , Φ i + ∇ • F(W h ), Φ i = -θ τ DDES (W h ), Φ i -(1 -θ) τ LES (W h ), Φ i
where W h denotes the hybrid variables and θ is the blending function defined as in the previous subsection.

NUMERICAL DISCRETIZATION

Numerical scheme

The governing equations are discretized in space using a mixed finite-volume/finite-element method applied to unstructured tetrahedral grids. The adopted scheme is vertex-centered, i.e. all degrees of freedom are located at the vertices. The diffusive terms are discretized using P1 Galerkin finite-elements on the tetrahedra, whereas finite-volumes are used for the convective terms. The numerical approximation of the convective fluxes at the interface of neighboring cells is based on the Roe scheme [START_REF] Roe | Approximate Riemann solvers parameters vectors and difference schemes[END_REF] with low-Mach preconditioning [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF]. To obtain second-order accuracy in space, the Monotone Upwind Scheme for Conservation Laws reconstruction method (MUSCL) [START_REF] Van Leer | Towards the ultimate conservative scheme. IV: a new approach to numerical convection[END_REF] is used, in which the Roe flux is expressed as a function of reconstructed values of W at each side of the interface between two cells. A particular attention has been paid to the dissipative properties of the resulting scheme, since this is a key point for its successful use in LES, and therefore in simulations performed with a hybrid turbulence model. The numerical dissipation provided by the scheme used in the present work is made of sixth-order space derivatives [START_REF] Camarri | A low diffusion MUSCL scheme for LES on unstructured grids[END_REF] and thus is concentrated on a narrow-band of the highest resolved frequencies. This is expected to limit undesirable damping of the large scales by numerical dissipation. Moreover, a parameter γ directly controls the amount of introduced numerical viscosity and can be explicitly tuned in order to reduce it to the minimal amount needed to stabilize the simulation. Time advancing is carried out through an implicit linearized method, based on a second-order accurate backward difference scheme and on a firstorder approximation of the Jacobian matrix [START_REF] Martin | A secondorder defect correction scheme for unsteady problems[END_REF]. The resulting numerical discretization is second-order accurate both in time and space. It should be noted that the spatial discretization used in this work leads to a superconvergent approximation, i.e. the accuracy can be well above second-order for some Cartesian meshes. One can also add that the objective is not a higher-order convergence but a strong reduction of dissipation and a certain reduction of the dispersion in the general case of a non-Cartesian (but not too irregular) mesh.

Mesh adaptation for rotating machines

Our numerical model has been extended to rotorstator simulation with a Chimera technique.

The Chimera method aims at solving partial differential equations by decomposition into subdomains with overlap in order to avoid having to use a global mesh. This method allows the communication between the computational subdomains thanks to the overlapping of the subdomains. In our case we consider a decomposition in two domains as shown in Fig. 4, a fixed domain in red and a rotating domain in blue. A mesh adaptation loop based on a metric depending of the Mach number has been combined with the Chimera algorithm. The Transient Fixed Point mesh adaptation algorithm [START_REF] Alauzet | 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations[END_REF][START_REF] Alauzet | Multi-Scale Anisotropic Mesh Adaptation for Time-Dependent Problems[END_REF][START_REF] Guégan | An L ∞ -L p mesh adaptive method for computing unsteady bi-fluid flows[END_REF] is applied to each domain from a metric field evaluated on the whole domain.

APPLICATIONS

Flow past a cylinder

The predictions of the flow around a circular cylinder are presented. Two Reynolds numbers, 10 6 and 2 × 10 6 , based on the cylinder diameter, D, and on the freestream velocity, are considered. Only a few numerical investigations have been performed for Reynolds numbers higher than 10 6 . This interval is inside the supercritical regime which appears at Reynolds number higher than 2 × 10 5 and for which the separation becomes turbulent [START_REF] Sumer | Hydrodynamics Around Cylindrical Structures[END_REF]. The computational domain is such that -15 ≤ (x, y)/D ≤ 15, and -1 ≤ z/D ≤ 1 where x, y and z denote the streamwise, transverse and spanwise directions respectively, the cylinder axis being located at x = y = z = 0. The mesh contains 4.8 millions nodes. The Steger-Warming conditions [START_REF] Steger | Flux vector splitting for the inviscid gas dynamic equations with applications to finite difference methods[END_REF] are imposed at the inflow and outflow. The flow is assumed to be periodic in the spanwise direction. The inflow Mach number is set to 0.1 so that the compressibility effects can be considered as negligible. The freestream turbulence intensity is below 1%. In order to control the computational costs, a wall law is applied in the close vicinity of the wall. For accuracy purpose, the Reichardt analytical law [START_REF] Hinze | Turbulence[END_REF], which gives a smooth matching between linear, buffer and logarithmic regions, is chosen. As the y + normalized distance is generally subject to large variations in complex flows, the wall law is combined with low Reynolds modeling which locally damps the fully turbulent model in regions in which the wall law does not cover the buffer zone.

Reynolds number 10 6

The behavior of the hybrid models presented in Section 2 are first investigated in terms of flow bulk coefficients. From Table 1, it can be noted that the predictions of these coefficients are globally in good agreement with the experimental data and the numerical results in the literature. The lift time fluctuations are nevertheless better predicted by the RANS/DVMS model compared to the other two hybrid models. This good behavior is confirmed by the correct prediction of the distribution of the mean pressure coefficient, see Fig. 5. On the other hand, it can be observed from Table 1 and Fig. 5 that Smagorinsky and WALE SGS models give globally comparable results. It is also worth noting that the RANS model predicts rather well the bulk coefficients. Fig. 6 shows instantaneous isocontours of the vorticity magnitude for each of the turbulence models used for this benchmark. As one might expect, the unsteady RANS model is more dissipative and captures less flow detail than the hybrid models, with in particular a much more damped and regular wake. The main outputs obtained by the RANS, DDES, DDES/DVMS and RANS/DVMS models are summed up in Table 2. Regarding the two hybrid models using the DVMS approach, only the results obtained with the Smagorinsky SGS model are shown, the WALE SGS model giving very similar results. It can be noticed that all models overpredict the separation angle, and that the fluctuations of the lift coefficient are overpredicted by the RANS and DDES/DVMS model. The RANS/DVMS approach gives on the whole the most satisfactory results. In Fig. 7, which shows the distribution of the mean pressure coefficient, the numerical results obtained with both hybrid models RANS/DVMS and DDES/DVMS are in very good agreement with the experimental results, in particular RANS/DVMS.

Flow around a rotating cross

We present a preliminary computation obtained with the combination of the mesh adaptation, the Chimera approximation and the DDES model. They are used for simulating the mixing process obtained with the rotation of a cross in a cylindric box. The Reynolds number is 1.8 million, based on the thickness of the blades. The number of vertices of the initial mesh is 10.466. After two turns, and five mesh adaptations, the number of vertices of the adapted mesh is 1.05 million. In Fig. 8a, we present a view of this final adapted mesh. Mesh size normal to the cylinder is yet not small enough to allow the creation of strong non-2D features. This is verified from the examination of the Q-factor, Fig. 8b. Fig. 9 shows details of mesh and velocity magnitude. It seems from this first calculation that despite the hign Reynold number, the flow separates at blade angle is a rather stable mode.

CONCLUSION

Several hybrid strategies, based on the DDES, RANS and DVMS models, are first evaluated for the simulation of the supercritical flow around a circular cylinder. This benchmark is characterized by turbulent boundary-layer separation. The cylinder flow is first computed at Reynolds 10 6 . The predictions of the main flow parameters are in overall good agreement with experimental data and the results of simulations in the literature, especially the RANS/DVMS model. A second series of computations of the flow past a circular cylinder is carried out at Reynolds 2 × 10 6 . The predictions obtained by the hybrid models used to perform this benchmark are a little less accurate compared to the previous Reynolds number, more specifically for the mean separation angle which is overestimated. Again, the RANS/DVMS model is the hybrid model among those applied which gives the best results overall. It is also to be noted that the pressure distribution is well predicted by the hybrid models for both Reynolds numbers, especially by the RANS/DVMS model. Then, the DDES formulation is applied to the flow around a rotating cross inside a cylinder with a Reynolds number of 1.8 × 10 6 . Preliminary results using mesh adaptation with meshes of 1 million vertices show that separations are more stable than in the flows around the cylinder. RANS/DVMS computations are planned. The new tool is also being applied to other rotating devices like propellers and will produce mesh adapted results for these geometries in a near future.

Figure 1 .

 1 Figure 1. Building the VMS coarse level

Figure 2 .

 2 Figure 2. VMS principle.

  by replacing the D RANS k = ρε dissipation term in the k equation by D DDES k = ρ k 3/2

Figure 3 .

 3 Figure 3. Isocontours of the blending function for the circular cylinder benchmark (Re= 10 6 ):for θ = 1 (in red) the RANS model is activated, wherever 0 < θ < 1 additional resolved fluctuations are computed, and in the limit θ → 0 (in blue) the LES approach (DVMS in this work) is recovered.

  For computations we start by locating the boundary nodes of domain 1 in domain 2, and reciprocally we locate the boundary nodes of domain 2 in domain 1, then the aerodynamic values of each boundary node are determined by interpolation and each domain performs its calculation with the new interpolated values.

Figure 4 .

 4 Figure 4. Definition of rotor and stator from an initial mesh : in red, the stator, in blue the rotor, and in gray, the overlap region.

Figure 5 .

 5 Figure 5. Distribution of the mean pressure coefficient over the cylinder surface. Comparison between experimental data and numerical results at Reynods number 10 6 .

Figure 6 .

 6 Figure 6. Circular cylinder, Re= 10 6 : instantaneous isocontours of the vorticity magnitude. From top to bottom: RANS, DDES, DDES/DVMS and RANS/DVMS.

Figure 7 .

 7 Figure 7. Distribution of the mean pressure coefficient over the cylinder surface. Comparison between experimental data and numerical results at Reynods number 2 × 10 6 .

Figure 8 .

 8 Figure 8. DDES computation of a rotating mixing device with a mesh adaptative Chimera approach (horizontal cut). View of the mesh and of Q-factor (colored with velocity magnitude) after the vertical upper blade.

Figure 9 .

 9 Figure 9. DDES computation of a rotating mixing device with a mesh adaptative Chimera approach (horizontal cut). Partial view of the mesh and velocity magnitude.

Table 1 .

 1 Bulk coefficients of the flow around a circular cylinder at Reynolds number 10 6 . The subscripts S and W holds respectively for Smagorinsky and WALE SGS models.

		C d	C l	-C pb	θ
	Present simul.
	RANS	0.20	0.02	0.20	130
	DDES	0.20	0.04	0.22	138
	DDES/DVMS S	0.20	0.02	0.22	135
	DDES/DVMS W	0.20	0.02	0.26	132
	RANS/DVMS S	0.25	0.09	0.25	132
	RANS/DVMS W	0.26	0.11	0.22	134
	Other simul.
	RANS [26]	0.40	-	0.41	-
	LES [26]	0.31	-	0.32	-
	LES [27]	0.27	0.12	0.28	108
	Measurements
	Exp. [28]	-	0.24	0.33	-
	Exp. [29]	0.22	-	-	-
	Exp. [30]	0.25	0.32	-	-
	Exp. [31]	-	-	-	130
	Exp. [32]	0.2-0.4 0.1-0.15 0.2-0.34	-

Table 2 .

 2 Bulk coefficients of the flow around a circular cylinder at Reynolds number 2 × 10 6 . The subscript S holds for Smagorinsky SGS model.

		C d	C l	-C pb	θ
	Present simul.
	RANS	0.26 0.066	0.30	128
	DDES	0.28 0.038	0.27	132
	DDES/DVMS S	0.26 0.070	0.35	130
	RANS/DVMS S	0.24 0.030	0.30	132
	Other simul.
	LES/TBLE [33] 0.24 0.029	0.36	105
	Measurements
	Exp. [28]	0.26 0.033	0.40	105
	Exp. [29]	0.32 0.029	-	-

ACKNOWLEDGEMENTS

This work was supported by the ANR NORMA project, grant ANR-19-CE40-0020-01 of the French National Research Agency. The authors gratefully acknowledge GENCI for granted access to HPC resources through CINES (grants 2020-A0092A05067 and 2021-A0102A06386) and IDRIS (grant 2021-A0112A05067).