Pallavi Singh

Manuela Lopez Villavicencio

Steve Labrie
email: steve.labrie@fsaa.ulaval.ca

Vincent Perkins

Stéphanie Vignola

Marie-Hélène Lessard

Pier-Luc Plante

Jacques Corbeil

Eric Dugat-Bony

Michel Frenette

Phenotypic and Genetic Characterization of the Cheese Ripening Yeast Geotrichum candidum

Keywords: yeast, cheese, Geotrichum candidum, Galactomyces candidus, multilocus sequence typing, assimilation of carbon compounds, genome assembly

Genetic Characterization of Geotrichum candidum of 11 G. candidum/Galactomyces spp. genomes allowed us to optimize a genotyping MLST scheme and, combined with the assessment of their ability to grow under different conditions, provides a reliable tool to cluster and eventually improves the selection of G. candidum strains.

The yeast Geotrichum candidum (teleomorph Galactomyces candidus) is inoculated onto mold-and smear-ripened cheeses and plays several roles during cheese ripening. Its ability to metabolize proteins, lipids, and organic acids enables its growth on the cheese surface and promotes the development of organoleptic properties. Recent multilocus sequence typing (MLST) and phylogenetic analyses of G. candidum isolates revealed substantial genetic diversity, which may explain its strain-dependant technological capabilities. Here, we aimed to shed light on the phenotypic and genetic diversity among eight G. candidum and three Galactomyces spp. strains of environmental and dairy origin. Phenotypic tests such as carbon assimilation profiles, the ability to grow at 35 • C and morphological traits on agar plates allowed us to discriminate G. candidum from Galactomyces spp. The genomes of these isolates were sequenced and assembled; whole genome comparison clustered the G. candidum strains into three subgroups and provided a reliable reference for MLST scheme optimization. Using the whole genome sequence as a reference, we optimized an MLST scheme using six loci that were proposed in two previous MLST schemes. This new MLST scheme allowed us to identify 15 sequence types (STs) out of 41 strains and revealed three major complexes named GeoA, GeoB, and GeoC. The population structure of these 41 strains was evaluated with STRUCTURE and a NeighborNet analysis of the combined six loci, which revealed recombination events between and within the complexes. These results hint that the allele variation conferring the different STs arose from recombination events. Recombination occurred for the six housekeeping genes studied, but most likely occurred throughout the genome. These recombination events may have induced an adaptive divergence between the wild strains and the cheesemaking strains, as observed for other cheese ripening fungi. Further comparative genomic studies are needed to confirm this phenomenon in G. candidum. In conclusion, the draft assembly

INTRODUCTION

The yeast species Geotrichum candidum (teleomorph Galactomyces candidus) is commonly found in water, air, soil, cereals, ripened fruits, milk, and especially on the surface of mold-and smear-ripened cheeses like Camembert, Tilsit, and Pont-L'Évêque [START_REF] Marcellino | Diversity of Geotrichum candidum strains isolated from traditional cheesemaking fabrications in France[END_REF][START_REF] Boutrou | Interests in Geotrichum candidum for cheese technology[END_REF][START_REF] Thornton | Detection of the sour-rot pathogen Geotrichum candidum in tomato fruit and juice by using a highly specific monoclonal antibody-based ELISA[END_REF][START_REF] De Hoog | Galactomyces Redhead & Malloch (1977)[END_REF][START_REF] Eliskases-Lechner | Yeasts and molds | Geotrichum candidum[END_REF][START_REF] Desmasures | CHEESE | Mold-ripened varieties[END_REF]. G. candidum was long considered an undesirable microorganism because its presence in cheese was uncontrolled, but is now recognized for its contribution to cheese ripening, due to its capacity to consume lactic acid and generate alkaline products [START_REF] Marcellino | Diversity of Geotrichum candidum strains isolated from traditional cheesemaking fabrications in France[END_REF][START_REF] Spinnler | Surface mould-ripened cheeses[END_REF][START_REF] Boutrou | Interests in Geotrichum candidum for cheese technology[END_REF]. This yeast also contributes to the development of flavor in the cheese due to its proteolytic and lipolytic activities [START_REF] Boutrou | Contribution of Geotrichum candidum to the proteolysis of soft cheese[END_REF][START_REF] Gente | Levels of cystathionine γ lyase production by Geotrichum candidum in synthetic media and correlation with the presence of sulphur flavours in cheese[END_REF]. Additionally, G. candidum is a good competitor against unwanted microorganisms such as Mucor spp. and has an antagonistic action against pathogens such as Listeria monocytogenes [START_REF] Dieuleveux | Purification and characterization of anti-Listeria compounds produced by Geotrichum candidum[END_REF][START_REF] Boutrou | Interests in Geotrichum candidum for cheese technology[END_REF]. However, these various traits, especially those contributing to the development of cheese flavors and texture, are strain dependent. G. candidum species display a great diversity at the morphologic and metabolic levels [START_REF] Guéguen | Les levures et Geotrichum candidum[END_REF][START_REF] Jollivet | Production of volatile compounds in model milk and cheese media by eight strains of Geotrichum candidum Link[END_REF][START_REF] Jacobsen | Specific populations of the yeast Geotrichum candidum revealed by molecular typing[END_REF][START_REF] Molimard | Effect of the association of surface flora on the sensory properties of mould-ripened cheese[END_REF]. Their high level of polymorphism and widely variable phenotypic traits complicate the identification and characterization of G. candidum strains by traditional microbiological methods [START_REF] Prillinger | Phenotypic and genotypic identification of yeasts from cheese[END_REF][START_REF] Gente | Identification of Geotrichum candidum at the species and strain level: proposal for a standardized protocol[END_REF]. Many indigenous strains of G. candidum have been isolated from the environment or from dairy products [START_REF] Marcellino | Diversity of Geotrichum candidum strains isolated from traditional cheesemaking fabrications in France[END_REF]Gente et al., 2002b;[START_REF] Lavoie | Characterization of the fungal microflora in raw milk and specialty cheeses of the province of Quebec[END_REF]. However, the lack of genomic and functional information has limited their potential utility.

Different G. candidum strains can be discerned based on chromosome sizes or sequences of their housekeeping genes (Gente et al., 2002a;[START_REF] Alper | Genetic diversity of dairy Geotrichum candidum strains revealed by multilocus sequence typing[END_REF][START_REF] Morel | Differential gene retention as an evolutionary mechanism to generate biodiversity and adaptation in yeasts[END_REF][START_REF] Jacobsen | Specific populations of the yeast Geotrichum candidum revealed by molecular typing[END_REF]. The ability to identify distinct phylogenetic groups within the G. candidum species lead to the development of two independent multilocus sequence typing (MLST) schemes to genotype isolates of G. candidum, but those were never compared using the same set of strains [MLST2013 [START_REF] Alper | Genetic diversity of dairy Geotrichum candidum strains revealed by multilocus sequence typing[END_REF]) and MLST2017 (Jacques et al., 2017)]. Moreover, only one genome of the G. candidum species has been partially sequenced and subjected to further gene prediction analysis [START_REF] Morel | Differential gene retention as an evolutionary mechanism to generate biodiversity and adaptation in yeasts[END_REF].

The aim of this study was to combine phenotypic assays and whole genome sequencing to elucidate the phenotypic and genetic diversity of eight G. candidum strains that belong to distantly related groups from a previous study using MLST [START_REF] Alper | Genetic diversity of dairy Geotrichum candidum strains revealed by multilocus sequence typing[END_REF]; Table 1). We also used the newly sequenced genomes to optimize a MLST scheme based on the two MLST schemes available in the literature [START_REF] Alper | Genetic diversity of dairy Geotrichum candidum strains revealed by multilocus sequence typing[END_REF][START_REF] Jacobsen | Specific populations of the yeast Geotrichum candidum revealed by molecular typing[END_REF]. The phenotypic and genomic results obtained in this study increase the scientific knowledge on the ripening yeast G. candidum and, in the near future, could be used to improve strain selection for cheese ripening processes.

MATERIALS AND METHODS

Biological Material, Culture Conditions, and Genomic DNA Extraction

Eight G. candidum (identified in bold in Table 1) and three Galactomyces spp. strains were used in this study (Table 2). The eight G. candidum strains were selected as representatives of each branch in the phylogenetic tree obtained with MLST2013 [START_REF] Alper | Genetic diversity of dairy Geotrichum candidum strains revealed by multilocus sequence typing[END_REF], while the Galactomyces spp. strains were chosen to form an outgroup of other Galactomyces species. Considering that the genus Galactomyces is not well known, the choice of three strains of different species as outgroups allowed to draw a wider portrait of this genus. All strains were grown on YEG (Yeast Extract Glucose) agar plates [10 g•L -1 of yeast extract (Thermo Fischer Scientific), 10 g•L -1 of D-glucose (EMD Chemicals) and 15 g•L -1 of Bacto agar (BD Diagnostics)]. The strains were inoculated directly from 15% glycerol (v/v) stock cultures stored at -80 • C. The plates were incubated in the dark for 4 days at 25 • C. Isolated colonies were re-streaked on YEG agar plates or in YEG Broth and incubated for 7-9 days at 25 • C in the dark. An additional 33 G. candidum strains isolated from dairy environment was used to validate the newly designed MLST scheme (Table 1).

The mycelia on YEG Broth were harvested and grounded into fine powder in liquid nitrogen using a CryoMill apparatus (Retsch, Germany). Cryogenic grinding was performed at -196 • C with an automatic pre-cooling step at a 5 Hz frequency followed by a 2 min grinding step at a 25 Hz frequency. Genomic DNA was extracted from 30 mg of grounded mycelium using Purelink RNA/DNA viral mini kit (Invitrogen) with the following modifications. Grounded mycelium was homogenized in 200 µl of a 0.9% NaCl solution, prior to the addition of the proteinase K and the lysis buffer. DNA concentration and quality were measured using a NanoDrop ND-1000 spectrophotometer The arthrospores of eight G. candidum and three Galactomyces spp. strains were isolated from the mycelium with a sterile cotton swab using 0.5% (v/v) Tween 80, suspended in YNB starvation broth [6.7 g•L -1 of Bacto Yeast Nitrogen Base without amino acids (BD Diagnostics) and 100 g•L -1 of D-glucose (EMD Chemicals)], and incubated for 48 h at 25 • C in the dark [START_REF] Kurtzman | Methods for isolation, phenotypic characterization and maintenance of yeasts[END_REF]. Microscopic images were taken, and spore concentration was determined using a hemocytometer and optical microscopy (BX61, Olympus). Thereafter, the solution was diluted with sterile carbon-free YNB solution to obtain a concentration of 1•10 6 spores•ml -1 . In order to compare the morphology between strains, isolated colonies were re-streaked in three points on YEG agar plates and photographed after 7 days of incubation at 25 • C using a standardized photography setup [START_REF] Pitt | Fungi and Food Spoilage[END_REF]. Their ability to grow at 35 • C was evaluated by inoculating G. candidum/Galactomyces spp. strains on GYP agar plates [40 g•L -1 of D-glucose (EMD Chemicals), 5 g•L -1 of yeast extract (Thermo Fischer Scientific), 5 g•L -1 of Bacto peptone (BD Diagnostics) and 20 g•L -1 of Bacto agar (BD Diagnostics)] for 4 days at 35 • C in sealed plastic bags (containing O 2 ) immersed in a water bath to ensure a constant temperature [START_REF] Kurtzman | Methods for isolation, phenotypic characterization and maintenance of yeasts[END_REF]. Plates were inoculated in a single point with 1 µl from the starved and diluted yeast cell suspension (1•10 6 spores•ml -1 ). Growth was assessed by the presence or the absence of a colony. Preliminary screening of carbon compound assimilation abilities was performed on the same 11 G. candidum/Galactomyces spp. strains using YT MicroPlates (BIOLOG) according to the manufacturer's protocol, with the following modifications. All wells of the microplate were filled with 100 µl of the starved and diluted yeast cell suspension (1•10 6 spores.ml -1 ). Plates were incubated in the dark for 10 days at 25 • C in unsealed plastic bags to avoid drying out of the wells. Yarrowia lipolytica LMA-23 was used as a positive control, while the negative control consisted of sterile distilled water. Yeast growth was determined after 3 and 10 days by optic density (OD) measurement with a BioTek microplate reader Synergy HI and Gen5 TM software v2.07. A scan area measurement of the OD was done at a wavelength of 590 nm. Carbon compounds, namely adonitol (Acros Organics), D-arabitol (Acros Organics), D-(-)-amygdalin (Acros Organics), bromosuccinic acid (Sigma-Aldrich), α-ketoglutaric acid (Sigma-Aldrich), L-(-)-malic acid (Acros Organics), L-glutamic acid (MP Biomedicals), for which variations in assimilation profiles were observed between strains (data not shown) were selected for further tests. Additional carbon sources, namely lactose (Laboratoire MAT), D-glucono-1,5-lactone (Sigma-Aldrich), sodium citrate (BDH), sodium DL-lactate (Sigma-Aldrich) were also selected because their assimilation allows the formal identification of Geotrichum and Galactomyces strains at the species level and because of their relevance in the dairy environment [START_REF] De Hoog | Galactomyces Redhead & Malloch (1977)[END_REF][START_REF] Hill | Dairy products: cheese and yogurt[END_REF]. Each carbon source was added at 0.5% concentration (w/v) in 4.9 ml of YNB culture media. For acidic or basic compounds, the pH of the solutions was adjusted between 5.2 and 5.6, with HCl (1N) or NaOH (1N). Spores of all G. candidum and Galactomyces spp. strains (Tables 1,2) were then inoculated as described for the YT MicroPlates protocol, in triplicate, using the standard method of Wickerham and Burton and incubated for 21 days at 25 • C in the dark [START_REF] Wickerham | Carbon assimilation tests for the classification of yeasts[END_REF][START_REF] Kurtzman | Methods for isolation, phenotypic characterization and maintenance of yeasts[END_REF]. A 15 s vortex agitation was done on each test tube at days 7 and 21. Growth evaluation was determined by eyes using the descriptors established by [START_REF] Kurtzman | Methods for isolation, phenotypic characterization and maintenance of yeasts[END_REF].

Library Preparation, DNA Sequencing and Genome Assembly

Sequencing library construction was performed from 50 ng of total genomic DNA using the Illumina Nextera DNA library preparation kit (Illumina, San Diego, CA, United States) and Nextera index kit (Illumina), according to the manufacturer's instructions. Purified libraries were quantified using PicoGreen (Promega, Madison, WI, United States), diluted at 2 µM, multiplexed and sequenced on the Illumina HiSeq 2000 platform (101 bp paired-end reads). For strain LMA-244, the Illumina MiSeq platform (300 bp paired-end reads) was also used. FASTQ file generation and demultiplexing were performed using bcl2fastq v1.8.4 (Illumina, San Diego, CA, United States). Both single-and paired-end reads generated were quality filtered using FastQc v0.11.4 and trimmed with Trimmomatic v0.35 in order to cut adapters or indexes, to remove low quality bases (quality score <10) from the end of the reads, and to discard reads below 21 nt [START_REF] Andrews | FastQC: A Quality Control Tool for High Throughput Sequence Data[END_REF][START_REF] Bolger | Trimmomatic: a flexible trimmer for Illumina sequence data[END_REF]. Reads that remained unpaired were retained as single-end reads. Genome assembly was performed using SPAdes v3.11.1 [START_REF] Bankevich | SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing[END_REF]. The genome of strain CLIB 918 as reference was used for the assembly of the eight G. candidum genomes [START_REF] Morel | Differential gene retention as an evolutionary mechanism to generate biodiversity and adaptation in yeasts[END_REF] and de novo sequencing and assembly was performed for the Galactomyces strains. Based on the optimization of the total number of contigs, the number of contigs ≥1,000 bp and the N50 value, the k-mer option -k 21, 33, 55, 77 were used for all strains.

Scaffolds were filtered using the khmer software with a length cut-off of 1,000 bp [START_REF] Crusoe | The khmer software package: enabling efficient nucleotide sequence analysis[END_REF]. Genome assembly statistics and completeness were assessed using QUAST web interface and BUSCO v3 [START_REF] Gurevich | QUAST: quality assessment tool for genome assemblies[END_REF][START_REF] Waterhouse | BUSCO applications from quality assessments to gene prediction and phylogenomics[END_REF]. BUSCO measures the completeness of a genome assembly by comparing predicted genes in the assembly to a dataset of near-universal single-copy orthologs selected from OrthoDB v9 [START_REF] Waterhouse | BUSCO applications from quality assessments to gene prediction and phylogenomics[END_REF].

Genomic Content Comparison

A measure of the similarity between all the strains was generated by comparing the k-mer (31 nt) content (presence or absence) of the whole assembled genomes (scaffolds) using Ray Surveyor software, which is a functionality of Ray v3.0.0 [START_REF] Deraspe | Phenetic comparison of prokaryotic genomes using k-mers[END_REF]. Considering the set of all 31 bp k-mers of an assembled genome i, A i = (k 1 , k 2 , . . ., k lA ), and the set of all 31 bp k-mers of an assembled genome j, A j = (k 1 , k 2 ,. . ., k lB ), the set of all 31 bp k-mers shared between the assembled genomes i and j is defined as K i,j = (A i ∩ A j ). In order to have comparable values because of the varying genome sizes of the G. candidum and Galactomyces spp. strains, the normalized count of all 31 bp k-mers shared between the assembled genomes i and j, defined as

K i,j = |k i,j | √ |k ii |×|k jj |
, was rather used. Hence, a similarity heatmap was generated with gplots packages in the R environment using the values of the normalized set K i,j of shared 31 bp k-mers for all pairwise strain combinations [START_REF] Warnes | gplots: Various R Programming Tools for Plotting Data[END_REF]. These normalized set of values were also transformed in Euclidean distance metrics with the distance matrix computation function in the R environment and used to generate a UPGMA dendrogram with MEGA v6 [START_REF] Tamura | MEGA6: molecular evolutionary genetics analysis version 6.0[END_REF]R Development Core Team, 2015).

MLST Loci Selection

Two MLST schemes were developed and published independently for G. candidum isolates typing. MLST2013 targets six housekeeping genes coding for an alanyl-tRNA synthetase (ALA1), a pyruvate kinase (CDC19), an acetyl-CoA acetyltransferase (ERG10), a glutaminyl-tRNA synthase (GLN4), a phosphoglucoisomerase (PGI1), and a phosphoglucomutase (PGM2) [START_REF] Alper | Genetic diversity of dairy Geotrichum candidum strains revealed by multilocus sequence typing[END_REF]. MLST2017 targets five loci coding for a nucleoporin (NUP116), a dihydroorotate dehydrogenase (URA1), an oritidine-5 -phosphate decarboxylase (URA3), a yeast subtilisin-like protease III (SAPT4), and a phospholipase B (PLB3) [START_REF] Jacobsen | Specific populations of the yeast Geotrichum candidum revealed by molecular typing[END_REF].

All 11 loci targeted in MLST2013 and MLST2017 were amplified for the eight isolates of G. candidum (identified in bold in Table 1), using the primers available in Table 3 and the following parameters: initial denaturation step of 2 min at 94 • C, followed by 35 cycles (denaturation 30 s at 94 • C, annealing for 30 s at 56 • C and extension of 1 min at 72 • C) and a final extension of 3 min at 72 • C. PCR products were visualized by electrophoresis on 0.8% agarose gel, stained with Gel Red, and illuminated under UV light. For all primer sets, a single band was obtained. PCR fragments were sequenced using an ABI 3730×l Data Analyzer (Thermo Fisher Scientific, Waltham, MA, United States) at the Plate-forme de séquençage et de génotypage des génomes (Centre de recherche du CHU de Québec, Université Laval, Québec, Canada).

The dendrogram obtained with Ray Surveyor software was used as a reference to compare the phylogenetics trees constructed with the MLST2013 and MLST2017 schemes. Once the optimum consensus MLST scheme was designed, all six target loci (ALA1, CDC19, SAPT4, GLN4, PGI1, and PGM2) (Table 3) were amplified, using the same PCR conditions, on an additional set of 33 G. candidum strains (Table 1).

MLST Data Treatment and Bioinformatics Analysis

The Staden Package was used for alignment, edition and construction of consensus sequences based on the ABI sequence chromatograms [START_REF] Staden | The Staden sequence analysis package[END_REF]. Consensus sequences were entered in Unipro UGENE v1.31.1 software and multiple sequence alignments were performed using MUSCLE [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF][START_REF] Okonechnikov | Unipro UGENE: a unified bioinformatics toolkit[END_REF]. Sequences were trimmed to be in frame and encode the same number of amino acids. Evolutionary histories were inferred using the PhyML Maximum-likelihood method with bootstrap test of 1,000 replicates [START_REF] Felsenstein | Confidence limits on phylogenies: an approach using the bootstrap[END_REF].

For the new consensus MLST scheme (MLST2019), the number of polymorphic sites, the dN/dS ratio (dN is the number of non-synonymous substitutions per non-synonymous site and dS is the number of synonymous substitutions per synonymous site), and the average number of nucleotide differences between populations were determined using DnaSP software v6 [START_REF] Rozas | DnaSP 6: DNA sequence polymorphism analysis of large data sets[END_REF]. The allelic profile of each locus was determined based on nucleotide polymorphisms. The allelic profile of each locus was then assigned a unique number in order of discovery (e.g., Ala1-1, Ala1-2. . .). Afterward, each unique combination of six loci numbers (1-1-1-1-1-1) was used to determine the sequence type (ST; e.g., ST1) of each strain [START_REF] Jolley | Using MLST to study bacterial variation: prospects in the genomic era[END_REF].

Phylogenetic Analysis

The MLST2019 data were used to predict the population structure of G. candidum by using the STRUCTURE v2.3.4 software [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF]. Analyses were performed with a length of Burnin period of 80,000 and 80,000 rounds of calculation after the Burnin period. The number of assumed populations was determined by running simulations with K = 2 up to K = 10 with five iterations for each assumed population [START_REF] Porras-Hurtado | An overview of STRUCTURE: applications, parameter settings, and supporting software[END_REF], and was best estimated when the log (probability of data) ceased to increase rapidly. The population structure was then calculated with the optimal ancestral subpopulation (K = 5) with default parameters, including 30 iterations of the analysis for K = 5. Recombination analysis was performed for the whole data set, within and between the Geo complexes using the NeighborNet method and the phi-test for recombination implemented in the software SplitsTree v4 [START_REF] Bryant | Neighbor-net: an agglomerative method for the construction of phylogenetic networks[END_REF][START_REF] Huson | Application of phylogenetic networks in evolutionary studies[END_REF]. Additionally, the linkage disequilibrium between and within the loci was verified with the Index of Association (I A ) and the Standardized Index of Association (I A S ) using the LInkage ANalysis (LIAN) software v3.7 [START_REF] Smith | How clonal are bacteria?[END_REF][START_REF] Haubold | LIAN 3.0: detecting linkage disequilibrium in multilocus data[END_REF].

Data Availability

Raw sequence data for the Whole Genome Shotgun projects were deposited at GenBank under the Bioprojects accession numbers PRJNA482576 for LMA-40, PRJNA482605 for LMA-70, PRJNA482610 for LMA-77, PRJNA482613 for LMA-244, PRJNA482616 for LMA-317, PRJNA490507 for LMA-563, PRJNA482619 for LMA-1028, PRJNA490528 for LMA-1146, PRJNA486748 for LMA-1147, PRJNA486749 for LMA-1148, and PRJNA486756 for LMA-1150.

MLST partial gene sequences were deposited at GenBank under the accession numbers from MH745581 to MH745588 for ALA1, from MH745589 to MH745596 for CDC19, from MH745597 to MH745604 for ERG10, from MH745605 to MH745612 for GLN4, from MH745613 to MH745620 for NUP116, from MH745621 to MH745628 for PGI1, from

RESULTS

Colony Morphology

Macroscopic and microscopic observations on YEG culture media revealed three different morphotypes among the 11 G. candidum and Galactomyces spp. strains (i.e., yeast-like, intermediate or mold-like morphotypes) (Figure 1), as previously described [START_REF] Boutrou | Interests in Geotrichum candidum for cheese technology[END_REF]. Strains LMA-77, 563, 1028, and 1148 were characterized by mold-like, white, filamentous, resistant-to-the-touch and not-greasy colonies, with a predominance of vegetative hyphae, whereas LMA-1146 and 1150 were characterized by yeast-like, cream-colored, not resistant-to-the-touch and greasy colonies, with abundant production of arthrospores. Strains LMA-40, 70, 244, 317, and 1147 showed an intermediate morphology, characterized by white-felted, slightly resistant-to-the-touch and not-greasy colonies, with abundant production of arthrospores.

Carbon Assimilation Profile

These 11 G. candidum and Galactomyces spp. strains were previously identified by rDNA sequencing [START_REF] Alper | Ribosomal DNA polymorphisms in the yeast Geotrichum candidum[END_REF]

Sequencing, Genome Assembly and Genomic Content Comparison

The statistics for the eight G. candidum genome assemblies (Table 1 in bold) were highly variable, except for the assembly size, the GC content and the BUSCO completeness (Table 5).

Genome sizes varied from 23.1 Mb (LMA-40) to 24.2 Mb (LMA-70), with an average GC content of 42%. Similarly, the statistics for the three Galactomyces spp. de novo assemblies were highly variable depending on the species, except for the BUSCO completeness. The genome length varied between 22 Mb (Galactomyces geotrichum LMA-1147) and 26 Mb (Galactomyces citri-aurantii LMA-1150), and the GC content between 37% (Galactomyces reessii LMA-1148) and 40% (G. geotrichum LMA-1147). Although different strains and/or species were assembled, for all partial genomes, the BUSCO completeness percentage was close to 90% which suggests a good quality draft genome for all isolates sequenced. The similarity among the partially assembled genomes of the G. candidum and Galactomyces spp. strains was assessed by comparing their k-mer content, using the software implementation Ray Surveyor [START_REF] Deraspe | Phenetic comparison of prokaryotic genomes using k-mers[END_REF]. The dendrogram generated from the K i,j Euclidean distances and the heatmap representation of the normalized similarity values (K i,j) uncovered separate groups within the strains (Figure 2, left and center). First, the G. candidum strains (LMA-40, LMA-70, LMA-244, LMA-317, LMA-563, LMA-1028, and LMA-1146) were clearly separated from the other Galactomyces (strains LMA-1147, LMA-1148, and LMA-1150). These two groups shared 3% or less of their k-mer content. Even though G. candidum strains were closely related, they could be separated into three major groups based on their k-mer content (Figure 2).

Multilocus Sequence Typing Schemes for G. candidum

We observed similarities between both previous MLST schemes, according to the clustering of the strains (Supplementary Figures S1-S3). LMA-40, -77, -563, -1028, and-1146 were grouped together whereas LMA-70 and -244 were grouped together at the opposite side of the tree. On the other hand, LMA-317 clustered differently in the MLST2017 scheme compared to MLST2013. In MLST2017, LMA-317 clustered in a third group of isolates located between the other complexes. Also, for LMA-40, -77, -563, -1028, and-1146, clustering between strains of the same clade was different in the two schemes. This observation led us to evaluate the resolution level of both schemes (Supplementary Figures S1,S2) by constructing a maximum-likelihood phylogenetic tree for each scheme. We compared both trees to the dendrogram that we constructed using the distance matrix based on whole genome content (Figure 2 left panel). Both MLSTs showed similar clustering to the one using whole genomic content, but none had the same precision as the whole genome comparison. Following these observations, we selected six (6) target loci from both MLST schemes to develop an improved genotyping method for G. candidum isolates (Supplementary Figure S3), named MLST2019. We selected the loci ALA1, CDC19, SAPT4, GLN4, PGI1, and PGM2 because they clustered the isolates the same way as the shared k-mer comparison using the whole genomic content. We then genotyped 33 additional G. candidum strains using MLST2019.

We applied the MLST2019 scheme to a total of 41 strains of G. candidum (for which 40/41 were the same as MLST2013) and observed 15 STs distributed into three complexes, namely GeoA, GeoB, and GeoC (Figure 3). GeoA complex, supported 1). The GeoC complex was composed of strains isolated from various environments and were not clustered according to their isolation source. These strains had either a filamentous (LMA-21) or an intermediate (LMA-70, -244, and -1035) morphotype on YEG agar. LMA-1035 was a commercial ripening starter, and LMA-21, -70, and -244 were respectively isolated from rotting carrots, smear cheese and milk (Table 1).

The combined sequences of the six loci used in the MLST2019 scheme generated a concatenated contig of 2,859 bp with 64 polymorphic sites. The number of polymorphic sites varied from 5 for CDC19 to 18 for GLN4 and consisted mostly of synonymous substitutions. The number of non-synonymous substitution sites varied between 1 and 2 for most loci and up to 4 sites for the PGI1 locus. The allelic profile for each locus was 2 for SAPT4, 4 for GLN4 and PGI1, 5 for ALA1, 6 for CDC19, and 8 for PGM2.

The proportion of variable sites varied between 1.27 and 3.23%.

For each locus, the number of alleles was always lower than the number of polymorphic sites, which resulted in a number of variable alleles per variable sites ranging from 0.13 to 1.20. The dN/dS ratio was always close to 0 and mostly lower than 0.2, except for PGI1 (Table 6).

Evaluation of the Population Structure of G. candidum

Using MLST2019 ST profiles in STRUCTURE software, we statistically estimated the number of ancestral genotypes (K)

of the yeast G. candidum. The STs were clustered in 5 subpopulations (K = 5), called lineages or clades (1, 2, 3, 4, and 5) (Figure 4 and Table 1). When compared to MLST2019 clusters (GeoA, B, and C), lineages 2, 3, and 5 form the GeoA complex while lineage 1 corresponds to GeoC. The lineage 4 corresponds mostly to GeoB except for strains LMA-1026 and LMA-244, which come from GeoA and GeoC, respectively (Figure 3 and Table 1). The subpopulations observed did not seem to cluster the strains either by their ability to assimilate carbon sources, their colonial morphology or their isolation source, although lineages 2, 3, and 5 do not contain environmental isolates (Table 1 and Figure 4). All MLST2019 loci and Geo complexes were investigated for genetic exchanges by a NeighborNet analysis, based on nucleotide sequences. The six targeted loci (ALA1, CDC19, SAPT4, GLN4, PGI1, and PGM2) showed a tree-like subpopulation structure, which suggests that the allele diversity is due to clonal descent within a given locus rather than recombination (Supplementary Figure S4). The NeighborNet analysis was then performed for GeoA, GeoB, and GeoC complexes separately, FIGURE 4 | Estimation of the population structure of 41 G. candidum strains based on their sequence type profile. For each of the strains, on the abscissa, the estimated contribution, on the ordinate, to its genotype from each of five hypothetical ancestral genotypes is represented by bar plot of different colors. The source of isolation of each strain is shown above the histogram: cheese strain (C), environmental strain (E), industrial strain (I), milk (M). Each color represents a lineage in the G. candidum population. Lineage 1 (60% red and 40% yellow), lineage 2 (green), lineage 3 (blue), lineage 4 (close to 60% yellow and 40% red), and lineage 5 (purple).

for all pairs of complexes GeoA-B, GeoB-C, and GeoA-C, and for all complexes together GeoA-B-C (Figure 5). Tree-like structures were observed within GeoA, and GeoC complexes, and between GeoA and GeoC. The GeoB individual complex and all other comparisons (GeoA-B, GeoB-C, and GeoA-B-C) showed parallelogram-shaped (network-like) structures, suggesting intergenic recombination (data not shown). The recombination hypothesis was verified doing a phi-test in SplitsTree (Table 7; [START_REF] Bruen | A simple and robust statistical test for detecting the presence of recombination[END_REF]. The recombination events were supported statistically within the GeoB complex (p = 0.03), between the GeoA and the GeoB complexes (p ≤ 0.001) and on the whole data set (p ≤ 0.001) while the absence of recombination was statistically supported for the GeoA (p = 1.0), between GeoA and GeoC (p = 0.2) and between the GeoB and GeoC (p = 0.3). Because only four strains formed the GeoC complex, not enough data was available to perform the phi-test.

Finally, we verified the linkage disequilibrium of all loci within and between the Geo complexes by using the value of the Index of Association (I A ) and the Standardized Index of Association (I A S ) calculated using LIAN software. Both indexes measure the extent of linkage equilibrium within a population by quantifying the amount of intergenic recombination and by detecting association among alleles at different loci [START_REF] Agapow | Indices of multilocus linkage disequilibrium[END_REF]. If the value of the indexes differs significantly from zero, it is a sign of linkage disequilibrium and of a clonal population, while a value of 0 is a sign of linkage equilibrium due to frequent recombination [START_REF] Xu | Multilocus sequence typing of Lactococcus lactis from naturally fermented milk foods in ethnic minority areas of China[END_REF]. The I A S values were close to 0 for GeoB, GeoA-B, GeoB-C, and GeoA-B-C complexes, and were significantly different from 0 for GeoA (0.4, p = 0.5), GeoC (-0.13, p = 1), and GeoA-C (0.37, p ≤ 0.0001) complexes, confirming what was observed in the NeighborNet analysis (Table 7).

DISCUSSION

Phenotypic Identification and Characterization

The strains used in this study were successfully associated to G. candidum or a Galactomyces species following their phenotypic characterization [START_REF] De Hoog | Galactomyces Redhead & Malloch (1977)[END_REF]. The 4), as already observed (Kurtzman et al., 2011a). In addition, the pairwise comparison of the shared k-mer content (31-mers) confirmed the identification and phylogeny of the four species studied (Figure 2).

The G. candidum strains displayed variability in the assimilation of a few carbon compounds (Figure 2, right panel) which could reflect subgroups with different functional traits. In contrast with previous reports, we found that LMA-40 can use adonitol and did not detect variable assimilation of DL-lactate and D-mannitol among the G. candidum strains [START_REF] Marcellino | Diversity of Geotrichum candidum strains isolated from traditional cheesemaking fabrications in France[END_REF][START_REF] Alper | Genetic diversity of dairy Geotrichum candidum strains revealed by multilocus sequence typing[END_REF]. In addition, using traditional microbiological tests, we observed for the first time the variable assimilation of malic acid and α-ketoglutaric acid, the weak assimilation of amygdalin, and the positive assimilation of bromosuccinic acid in G. candidum. Some of the carbon compounds tested in this study can be metabolized during cheese ripening by different microorganisms, including G. candidum strains. Lactate assimilation by G. candidum during the early phase of cheese ripening has been described several times and is positively associated with deacidification of the cheese surface, which is mandatory for the growth of acidsensitive bacteria [START_REF] Fox | Enzymology of cheese ripening[END_REF][START_REF] Marcellino | Diversity of Geotrichum candidum strains isolated from traditional cheesemaking fabrications in France[END_REF][START_REF] Gripon | CHEESE | Mould-ripened cheeses[END_REF][START_REF] Leclercq-Perlat | Controlled production of Camembert-type cheeses. Part I: microbiological and physicochemical evolutions[END_REF][START_REF] Mcsweeney | Biochemistry of cheese ripening[END_REF][START_REF] Lessard | Metatranscriptome analysis of fungal strains Penicillium camemberti and Geotrichum candidum reveal cheese matrix breakdown and potential development of sensory properties of ripened Camembert-type cheese[END_REF][START_REF] Arfi | Importance of curd-neutralising yeasts on the aromatic potential of Brevibacterium linens during cheese ripening[END_REF][START_REF] Boutrou | Interests in Geotrichum candidum for cheese technology[END_REF][START_REF] Castellote | Investigation of Geotrichum candidum gene expression during the ripening of Reblochon-type cheese by reverse transcription-quantitative PCR[END_REF][START_REF] Dugat-Bony | Overview of a surface-ripened cheese community functioning by meta-omics analyses[END_REF]. Citrate metabolism in lactic acid bacteria leads to the production of CO 2 and diacetyl, which contributes to the typical aspect and flavor of Gouda, Cottage, Quark, and Cheddar cheeses [START_REF] Fox | Enzymology of cheese ripening[END_REF][START_REF] Molimard | Review: Compounds involved in the flavour of surface mould-ripened cheeses: origins and properties[END_REF][START_REF] Mcsweeney | Biochemistry of cheese ripening[END_REF]. However, citrate utilization as the sole carbon source by G. candidum strains was found to be weak or negative (Figure 2, right panel), suggesting that metabolites produced by G. candidum through citrate catabolism would not be major contributors to the typical flavor of surface-ripened cheese varieties [START_REF] Marcellino | Diversity of Geotrichum candidum strains isolated from traditional cheesemaking fabrications in France[END_REF][START_REF] Leclercq-Perlat | Controlled production of Camembert-type cheeses. Part I: microbiological and physicochemical evolutions[END_REF][START_REF] Lessard | Metatranscriptome analysis of fungal strains Penicillium camemberti and Geotrichum candidum reveal cheese matrix breakdown and potential development of sensory properties of ripened Camembert-type cheese[END_REF]. Glutamic acid, which can be consumed by all G. candidum strains tested, is produced through the transamination of L-methionine. This enzymatic reaction, in the presence of an amino-acid acceptor, usually the α-ketoglutaric acid, contributes to ammonia and aroma production by G. candidum strains during cheese ripening [START_REF] Brennan | Bacterial surface-ripened cheeses[END_REF][START_REF] Boutrou | Interests in Geotrichum candidum for cheese technology[END_REF][START_REF] Cholet | Gene expression and biochemical analysis of cheese-ripening yeasts: focus on catabolism of L-methionine, lactate, and lactose[END_REF][START_REF] Lessard | Metatranscriptome analysis of fungal strains Penicillium camemberti and Geotrichum candidum reveal cheese matrix breakdown and potential development of sensory properties of ripened Camembert-type cheese[END_REF][START_REF] Mcsweeney | Biochemistry of cheese ripening[END_REF]. Finally, no relationship was observed between the genetic population structure and the phenotypic characteristics.

Genome Comparisons and MLST Analysis

Before this study, the only genome sequence of G. candidum available in public databases was for the strain CLIB 918, which was isolated from Pont-L'Evêque cheese [START_REF] Morel | Differential gene retention as an evolutionary mechanism to generate biodiversity and adaptation in yeasts[END_REF]. The draft genome assemblies for eight additional G. candidum strains presented in this study contained from 836 to 1,488 scaffolds (≥1,000 bp) with N50 values between 29,977 and 137,000 bp (Table 4). The high number of short scaffolds observed might arise from the short reads produced through the Illumina HiSeq technology and the presence of repeated regions, which are abundant in eukaryotic genomes [START_REF] Paszkiewicz | De novo assembly of short sequence reads[END_REF][START_REF] Yandell | A beginner's guide to eukaryotic genome annotation[END_REF][START_REF] Kuo | Fungal genomics: sequencing and annotation[END_REF]. Using the G. candidum CLIB 918 draft genome sequence as a reference allowed major improvements in the scaffolding step for the G. candidum strains. The assembly size, N50 value and length of the longest scaffold increased with CLIB 918 as a reference, in comparison to assemblies without CLIB 918 (Supplementary Table S1). Overall, the three Galactomyces spp. draft genome assemblies provided similar estimations of the genome sizes and GC contents to those previously proposed using the thermal denaturation method [START_REF] Smith | Genome comparisons in the yeastlike fungal genus Galactomyces Redhead & Malloch[END_REF]. These assemblies contain close to 3,000 scaffolds (≥1,000 bp) with N50 value around 13,000 bp. The absence of an available reference genome for the scaffolding of Galactomyces spp. strains could explain the more fragmented assemblies, compared to those of G. candidum strains.

Sequencing the whole genome of eight G. candidum and three Galactomyces spp. strains allowed us to use the pairwise comparison of their shared k-mer content as the reference clustering method for optimizing current MLST schemes. Two previous MLST schemes were compared for the genotyping of a common set of eight G. candidum isolates [START_REF] Alper | Genetic diversity of dairy Geotrichum candidum strains revealed by multilocus sequence typing[END_REF][START_REF] Jacobsen | Specific populations of the yeast Geotrichum candidum revealed by molecular typing[END_REF]. Genotyping with MLST2013 gave rise to two groups of isolates, as previously obtained, whereas a third group was obtained with MLST2017 [START_REF] Alper | Genetic diversity of dairy Geotrichum candidum strains revealed by multilocus sequence typing[END_REF][START_REF] Jacobsen | Specific populations of the yeast Geotrichum candidum revealed by molecular typing[END_REF]. Although a few subgroups differ between MLST2013 and MLST2017, both schemes clustered LMA-70 and LMA-244 together and apart from the other G. candidum strains.

Considering that neither MLST2013 nor MLST2017 allowed the same sub-grouping as the shared k-mer content, a selection of six loci was used to create an optimized MLST scheme (MLST2019) that would increase the discrimination power between the isolates. The MLST2019 scheme uses the loci ALA1, CDC19, SAPT4, GLN4, PGI1, and PGM2, which, once concatenated, generate a sequence of 2,859 bp. The evolution rate was assessed using the dN/dS ratio for each locus (between 0.000 and 0.616) and suggests that they are not subjected to positive selection, making them suitable gene targets for MLST [START_REF] Shabbir | Genetic analysis of Newcastle disease virus from Punjab, Pakistan[END_REF][START_REF] Yu | Multilocus sequence typing of Streptococcus thermophilus from naturally fermented dairy foods in China and Mongolia[END_REF].

The evidence for genetic diversity could indicate an adaptive divergence, a population differentiation under selection, between the wild strains and the cheese ripening strains of G. candidum, as observed for other fungi used in the agri-food industry such as Penicillium roqueforti [START_REF] Dumas | Independent domestication events in the blue-cheese fungus Penicillium roqueforti[END_REF], Penicillium camemberti [START_REF] Cheeseman | Multiple recent horizontal transfers of a large genomic region in cheese making fungi[END_REF], and Saccharomyces cerevisiae [START_REF] Legras | Adaptation of S. cerevisiae to fermented food environments reveals remarkable genome plasticity and the footprints of domestication[END_REF]. This hypothesis is supported by the observation that most of the non-dairy strains are clustered in GeoB and GeoC complexes, whereas GeoA complex contains all industrial cheesemaking strains studied here, except for LMA-1035, and other isolates from milk or cheese. When analyzing the population structure of the G. candidum isolates, we observed that the Geo complexes (Figures 3,5) were similar, but there was a discordance in comparison to the estimation of the ancestral subpopulations (lineages) (Figure 4). The five G. candidum lineages were clustered accordingly to MLST2019, except for linage 4 which contains all GeoB strains and includes LMA-244 and LMA-1026. According to the concatenated MLST loci and the whole genome sequences, these strains are part of GeoC and GeoA complexes, respectively (Figures 3,5). The fact that LMA-244 and -1026 share alleles with strains of the complex GeoB could hint at a speciation of environmental strains (GeoB) to form new complexes GeoA and GeoC [START_REF] Dumas | Independent domestication events in the blue-cheese fungus Penicillium roqueforti[END_REF]. Additionally, the calculation of I A and I A S indices, both close to 0, supports the hypothesis that the strains studied are in linkage equilibrium and that recombination did play a role in increasing the STs diversity, which was statistically demonstrated, except for GeoB-C (Table 7; [START_REF] Jolley | Sequence type analysis and recombinational tests (START)[END_REF][START_REF] Cangi | Recombination is a major driving force of genetic diversity in the anaplasmataceae Ehrlichia ruminantium[END_REF]. In this particular case, the low number of strains (LMA-21, -70, -244, and -1035) in the GeoC complex might explain the discordance between both analysis. The isolation and genome sequencing of more strains clustering in GeoC would be needed to further characterize the species G. candidum. Moreover, a lower diversity was observed in the GeoA complex than in GeoB and GeoC, as represented by the fewer number of branches and their smaller size (Figure 3). A similar phenomenon has also been observed for domesticated P. roqueforti strains, for which cheese complexes had fewer and smaller branches compared to the other clusters [START_REF] Dumas | Independent domestication events in the blue-cheese fungus Penicillium roqueforti[END_REF]. Furthermore, GeoC isolates were distant from the GeoA/B complexes on a phylogenetic level, and isolates from the environment (LMA-21) and the dairy industries (LMA-70, -244, and -1235) were obtained for the GeoC complex. These findings suggest that the GeoC complex could represent a new pool of G. candidum strains with an interesting biotechnological potential for the development of cheese products with typical properties.

CONCLUSION

In conclusion, the aim of this study was to combine traditional microbiological techniques with new whole genome shotgun sequencing to illuminate the phylogenetic relationships among G. candidum and Galactomyces spp. strains of environmental and dairy origin. Thereby, we report the draft assembly of 11 Geotrichum/Galactomyces spp. strains, increasing the amount of genomic data available for G. candidum and other Galactomyces species, for which no genome sequence assemblies were available. These strains were identified based on their ability to grow under different conditions and we propose a new MLST scheme to optimally genotype G. candidum isolates, because it better represents the whole genomic content of this species and allows a clear visualization of its population structure. Further studies are needed to test the adaptive divergence hypothesis and future comparative genomics studies will shed light on the differences between GeoC and GeoA/B. Also, the sequencing of additional genome, especially from non-dairy niche, are likely to improve the new MLST scheme and give a better understanding of the diversity within each complex. Phenotypic and genomic tools such as whole genome sequencing and the MLST scheme substantially improve our knowledge on this ripening yeast and could ultimately allow a better selection and control of G. candidum strains throughout the ripening of cheeses. comparative genomics analysis. M-HL performed DNA libraries for HiSeq sequencing. VP performed DNA extraction for MLST schemes comparison, MLST analysis, and genome assembly and comparative genomics analysis. P-LP and JC helped with the development of the bioinformatic scripts used for genome analysis.

FIGURE 1 |

 1 FIGURE 1 | Morphology of Geotrichum candidum and Galactomyces spp. strains on YEG culture media. (Left) Macroscopic observation of colonies. (Right) Microscopic observation of arthrospores. (A) Mold-like (example: LMA-1028). (B) Intermediate (example: LMA-244). (C) Yeast-like (example: LMA-1146).

FIGURE 2 |

 2 FIGURE 2 | Genomic similarity of three Galactomyces spp. and eight G. candidum strains represented by a heatmap obtained from the Euclidean distance between the normalized set of k-mers (31 bp fragments) shared between the partial genomes of Galactomyces spp. and G. candidum strains. Phenetic tree analyses were conducted in MEGA6. (Right panel) Ability of G. candidum/Galactomyces spp. strains to grow on carbon compounds (+ for positive assimilation, -for negative assimilation, and W for weak assimilation).

FIGURE 3 |

 3 FIGURE 3 | Phylogram of the MLST2019 scheme (ALA1, CDC19, SAPT4, GLN4, PGI1, and PGM2) constructed with the PhyML Maximum Likelihood method with 1,000 bootstraps (bootstrap values are shown next to the main nodes) on 41 isolates of Geotrichum candidum. Phylogenetic analyses were performed with Unipro UGENE v1.31.1.

  strains LMA-1147 (G. geotrichum), LMA-1148 (G. reessii), and LMA-1150 (G. citri-aurantii) had the typical carbon assimilation profiles and the ability to grow at 35 • C. Based on the same criteria, all other isolates except LMA-1028 were positively identified as G. candidum. It was also possible to distinguish G. candidum strains from the Galactomyces species, based on their weak or positive utilization of α-ketoglutaric acid and malic acid (Table

FIGURE 5 |

 5 FIGURE 5 | NeighborNet analysis of the combined six loci (ALA1, CDC19, SAPT4, GLN4, PGI1, and PGM2) obtained from 41 G. candidum isolates. Numbering in the figure corresponds to sequence type. Parallelogram formation indicate recombination events. NeighborNet analysis was performed with SplitsTree version 4.

TABLE 1 |

 1 Origins and allelic profiles of the 41 G. candidum strains analyzed.

	Strain	Geo	Lineage	ST			Allele				Origin a	References
		complex
					ALA1 CDC19 SAPT4 GLN4 PGI1 PGM2
	15	A	3	1	1	1	1	1	1	1	Mold-ripened cheese, Canada	Alper et al., 2011
	20	A	2	2	1	2	1	1	2	2	Mold-ripened cheese, Canada	Alper et al., 2011
	21	C	1	3	2	3	2	2	2	3	Clotted carrot, Japan	Guéguen and Jacquet, 1982
	28	A	5	4	1	1	1	1	1	4	Mold-ripened cheese, Canada	Alper et al., 2011
	34	A	2	2	1	2	1	1	2	2	Mold-ripened cheese, Spain	Alper et al., 2011
	37	A	2	2	1	2	1	1	2	2	Smear cheese, France	Alper et al., 2011
	38	A	2	2	1	2	1	1	2	2	Mold-ripened cheese, France	Alper et al., 2011
	39	A	5	4	1	1	1	1	1	4	Smear cheese, France	Alper et al., 2011
	40	A	5	4	1	1	1	1	1	4	Smear cheese, France	Dieuleveux et al., 1997
	48	B	4	5	1	1	2	3	2	4	Bioreactor contaminant	Alper et al., 2011
	70	C	1	3	2	3	2	2	2	3	Smear cheese, Canada	Alper et al., 2011
	73	A	5	4	1	1	1	1	1	4	Smear cheese, France	Gente et al., 2002a,b, 2006
	74	B	4	6	1	4	2	1	3	5	Grass, France	Gente et al., 2002a,b, 2006
	75	B	4	7	1	5	2	1	4	4	Corn silage, France	Gente et al., 2002a,b, 2006
	76	A	5	8	1	1	1	1	2	4	Milk, France	Gente et al., 2002b, 2006
	77	A	3	1	1	1	1	1	1	1	Milk, France	Gente et al., 2002b, 2006
	244	C	4	9	3	3	2	4	2	6	Milk, Canada	Lavoie et al., 2012
	317	B	4	10	4	4	2	3	2	7	Organic milk, Canada	Lavoie et al., 2012
	436	A	2	2	1	2	1	1	2	2	Industrial strain-company A	Alper et al., 2011
	562	A	5	4	1	1	1	1	1	4	Milk, MPP 1, Canada	Lavoie et al., 2012
	563	A	5	4	1	1	1	1	1	4	Milk, Canada	Lavoie et al., 2012
	645	A	5	4	1	1	1	1	1	4	Smear cheese, Canada	Lavoie et al., 2012
	655	A	5	4	1	1	1	1	1	4	Mold-ripened cheese, Canada	Lavoie et al., 2012
	664	B	4	11	4	4	2	1	2	4	Milk, MPP 1, Canada	Lavoie et al., 2012
	690	B	4	12	5	6	2	1	4	7	Milk, MPP 3, Canada	Lavoie et al., 2012
	1024	A	3	13	5	1	1	1	1	1	Industrial strain-company B	Alper et al., 2011
	1025	A	3	13	5	1	1	1	1	1	Industrial strain-company B	Alper et al., 2011
	1026	A	4	14	1	1	1	1	2	8	Industrial strain-company B	Alper et al., 2011
	1028	A	2	2	1	2	1	1	2	2	Industrial strain-company A	Alper et al., 2013
	1031	A	3	13	5	1	1	1	1	1	Industrial strain-company A	Alper et al., 2013
	1032	A	3	13	5	1	1	1	1	1	Industrial strain-company A	Alper et al., 2013
	1033	A	3	13	5	1	1	1	1	1	Industrial strain-company A	Alper et al., 2013
	1034	A	2	2	1	2	1	1	2	2	Industrial strain-company A	Alper et al., 2013
	1035	C	1	15	3	4	2	2	2	3	Industrial strain-company A	Alper et al., 2013
	1036	A	2	2	1	2	1	1	2	2	Industrial strain-company A	Alper et al., 2013
	1037	A	3	13	5	1	1	1	1	1	Industrial strain-company A	Alper et al., 2013
	1038	A	5	4	1	1	1	1	1	4	Industrial strain-company A	Alper et al., 2013
	1039	A	5	4	1	1	1	1	1	4	Industrial strain-company A	Alper et al., 2013
	1040	A	5	4	1	1	1	1	1	4	Industrial strain-company A	Alper et al., 2013
	1041	A	3	1	1	1	1	1	1	1	Industrial strain-company A	Alper et al., 2013
	1146	A	3	1	1	1	1	1	1	1	Artisanal sheep cheese, Slovakia CBS 11176
	a MPP: Milk Production Plant. Bold: Strains sequenced in this study.

TABLE 2 |

 2 Identification and origins of the Galactomyces spp. strains.

	Species	Strain	Origin	References
	Galactomyces geotrichum	1147	Soil extract, Puerto Rico	CBS 772.71
	Galactomyces reessii	1148	Cold water retting, Indonesia	CBS 179.60
	Galactomyces citri-aurantii	1150	Orange orchard soil extract, California	CBS 175.89

TABLE 3 |

 3 Primer sequences of 11 MLST target loci.

	Locus	Primer sequence 5 → 3	Locus size (bp)	References
	ALA1	Fwd-GCTCTTCGTGAGGTTCTTGG	424	Alper et al., 2013
		Rev-ACCTCGTAGGCATCAGTGCT
	CDC19	Fwd-CGCCAGTCAGAGAAGGAATA	355	Alper et al., 2013
		Rev-GTCGACCTGGTTCTTGACAC
	ERG10	Fwd-AACACAACATTTCCCGTGAG	652	Alper et al., 2013
		Rev-AGAGCTTAGCGTCAGCACTGA
	GLN4	Fwd-TGTTCTCAGAGGGTTTCCTG	558	Alper et al., 2013
		Rev-CCACATCTGAGGATTGTCGT
	PGI1	Fwd-ACCGCTGAGACTCTTCGCAA	472	Alper et al., 2013
		Rev-CTCCATGGAAAGCTGCTGGA
	PGM2	Fwd-GAACGGTGTCTACGGTCTTG	548	Alper et al., 2013
		Rev-TCAATGTACAGACGGATGGTC
	NUP116	For-ACCGCTACAACTGGATTTGG	425	Jacques et al., 2017
		Rev-GAGACCTGTTTGAGGGCTTG
	URA1	For-CAAGCCAATTGTGCTGAGAA	465	Jacques et al., 2017
		Rev-GGTGTCGTAGGGCAGTTGAT
	URA3	For-GCCAAAAAGACCAACCTGTG	470	Jacques et al., 2017
		Rev-CCTCATCCATACGGTTCTGC
	SAPT4	For-ATCATTAACACCCCGGCATA	501	Jacques et al., 2017
		Rev-GTGTCACCAAGCAGAGCAAA
	PLB3	For-AAGAATATCTGGGATCTTTC	393	Jacques et al., 2017
		Rev-TGAAGAAGAAGTACCAAGAA
	MH745629 to MH745636 for PGM2, from MH745637 to
	MH745644 for PLB3, from MH745645 to MH745685 for SAPT4,
	from MH745686 to MH745693 for URA1, and from MH745694
	to MH745701 for URA3.

TABLE 4 |

 4 Carbon assimilation profile for the G. candidum/Galactomyces spp. strains.

	Genus	Geotrichum candidum

+: Positive assimilation; -: Negative assimilation; W: Weak assimilation; N/A: Not available. or yeast-like morphotype and these were divided into two sub-complexes, according to their isolation source. Dairy strains LMA-317, 664, and 690 isolated from milk produced in Canada clustered separately from the non-dairy strains LMA-48, -74, and -75, isolated respectively from a bioreactor, grass and corn silage (Table

TABLE 5 |

 5 Assembly statistics for the genome of Geotrichum candidum and Galactomyces spp. strains.

	LMA	Mean	Assembly	Total no. of	Scaffold	Longest	GC content	BUSCO
	strain	coverage	size (Mb)	scaffolds	N50 (bp)	scaffold (bp)	(%)	Completeness (%)
	40	63×	23.10	1,268	32,639	167,475	41.8	90.3
	70	70×	24.17	1,450	44,909	195,322	41.4	90.2
	77	58×	23.32	1,488	29,977	158,856	41.7	90.5
	244	114×	23.20	836	137,000	210,021	41.5	90.3
	317	19×	23.36	1,370	35,529	209,353	41.8	90.3
	563	49×	23.30	1,427	31,326	124,845	41.7	89.9
	1028	34×	23.42	1,332	36,124	191,607	41.7	90.5
	1146	54×	23.35	1,363	32,892	152,162	41.7	90.3
	1147	64×	21.96	3,121	13,090	102,582	40.4	89.8
	1148	105×	22.56	2,775	13,420	64,496	36.6	89.5
	1150	55×	25.98	2,668	16,289	69,389	39	90.9

TABLE 6 |

 6 Loci characteristics targeted in the consensus MLST scheme used to genotype 41 isolates of G. candidum.

	Locus	Locus size	Number of	Number of	Proportion of	Number of alleles per	dN/dS
		(bp)	polymorphic sites*	alleles	variable sites (%)	variable site
	ALA1	424	8 (1)	5	1.89	0.63	0.000
	CDC19	355	5 (1)	6	1.41	1.20	0.171
	SAPT4	501	16 (2)	2	3.19	0.13	0.079
	GLN4	558	18 (2)	4	3.23	0.23	0.028
	PGI1	472	6 (4)	4	1.27	0.67	0.616
	PGM2	548	11 (1)	8	2.01	0.73	0.048
	*Numbers in parentheses are the number of non-synonymous substitution sites detected in the locus.

TABLE 7 |

 7 Recombination and linkage disequilibrium statistics.

	Geo group	Phi	I A	I A	S	p-value
	A	1.0	0.21	0.4	0.5
	B	0.03	-0.33	-0.07	0.9
	C	NA	-0.63	-0.13	1
	A-B	≤0.001	0.32	0.06	0.05
	B-C	0.3	-0.2	-0.04	0.8
	A-C	0.2	1.9	0.37	0.0001
	A-B-C	≤0.001	0.36	0.07	0.03

Frontiers in Microbiology | www.frontiersin.org

May 2020 | Volume 11 | Article 737

---------

ACKNOWLEDGMENTS

The realization of this work could not have been possible without the great help of Dr. Frédéric Raymond and the technical support of the Calcul Québec team, Félix-Antoine Fortin and Maxime Boissonneault. Computations were made on the supercomputer Cedar from the Simon Fraser University, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), Ministère de l'Économie, des Sciences et de l'Innovation du Québec (MESI) and le Fonds de recherche du Québec -Nature et technologies (FRQ-NT).

FUNDING

This work was supported by a Discovery grant from the Natural Sciences and Engineering Research Council of Canada (NSERC). Vincent Perkins was the recipient of the of a Levier Doctoral scholarship awarded from the INITIA Foundation. SV was the recipient of the Québec Graduate Scholarship (FRQ-NT) and the Canada Graduate scholarship (NSERC). Eric Dugat-Bony was supported by the AgreenSkills+ fellowship program which has received funding from the EU's Seventh Framework Programme under grant agreement N • FP7-609398 (AgreenSkills+ contract).

DATA AVAILABILITY STATEMENT

The datasets generated for this study can be found in the raw sequence data for the Whole Genome Shotgun projects were deposited at GenBank under the Bioprojects accession numbers PRJNA482576 for LMA-40, PRJNA482605 for LMA-70, PRJNA482610 for LMA-77, PRJNA482613 for LMA-244, PRJNA482616 for LMA-317, PRJNA490507 for LMA-563, PRJNA482619 for LMA-1028, PRJNA490528 for LMA-1146, PRJNA486748 for LMA-1147, PRJNA486749 for LMA-1148, and PRJNA486756 for LMA-1150. MLST partial gene sequences were deposited at GenBank under the accession numbers from MH745581 to MH745588 for ALA1, from MH745589 to MH745596 for CDC19, from MH745597 to MH745604 for ERG10, from MH745605 to MH745612 for GLN4, from MH745613 to MH745620 for NUP116, from MH745621 to MH745628 for PGI1, from MH745629 to MH745636 for PGM2, from MH745637 to MH745644 for PLB3, from MH745645 to MH745685 for SAPT4, from MH745686 to MH745693 for URA1, and from MH745694 to MH745701 for URA3.

AUTHOR CONTRIBUTIONS

VP, SV, M-HL, ED-B, MF, and SL were involved in planning the study and writing the manuscript. SV performed phenotypic characterization, DNA extraction, and genome assembly and

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2020.00737/full#supplementary-material

Conflict of Interest:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Perkins, Vignola, Lessard, Plante, Corbeil, Dugat-Bony, Frenette and Labrie. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.


Mario Verdugo-Marchese
email: mario.verdugo-marchese@chuv.ch.

MSc Stefano Coiro

MD, MSc Christine Selton-Suty

Masatake Kobayashi

Erwan Bozec

PhD Zohra Lamiral

Clément Venner

Faiez Zannad

PhD Patrick Rossignol

PhD Nicolas Girerd

MD Msc ; Faiez Zannad

PhD ; Olivier Huttin

Left ventricular myocardial deformation pattern, mechanical dispersion, and their relation with electrocardiogram markers in the large population-based STANISLAS cohort: insights into electromechanical coupling

Keywords: Myocardial deformation, Left ventricle, systolic function, Speckle-tracking echocardiography, Population study, Mechanical dispersion, Electromechanical coupling

à 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

Electrical conduction in the heart and cardiomyocyte depolarization culminates in mechanical contraction (excitation-contraction coupling) [START_REF] Pfeiffer | Biomechanics of cardiac electromechanical coupling and mechanoelectric feedback[END_REF]. Pathological alterations in magnitude and timing of myocardial activity result in alterations in electrocardiogram (ECG) recordings. These subtle electrical variations should logically be translated into mechanical variations.

Left ventricular ejection fraction (LVEF) has long been the main biomarker for sudden cardiac death (SCD) risk assessment and therapeutic decision-making [START_REF] Priori | ESC Guidelines for the management of patients with ventricular arrhythmias and the prevention of sudden cardiac death: The Task Force for the Management of Patients with Ventricular Arrhythmias and the Prevention of Sudden Cardiac Death of the European Society of Cardiology (ESC). Endorsed by: Association for European Paediatric and Congenital Cardiology (AEPC)[END_REF] although its accuracy has recently been questioned [START_REF] Seegers | ICD risk stratification studies -EU-CERT-ICD and the European perspective[END_REF][START_REF] Køber | Defibrillator Implantation in Patients with Nonischemic Systolic Heart Failure[END_REF]. Intraventricular conduction (QRS) delay and dispersion are linked to an increased risk of arrhythmia in a myriad of heart diseases [START_REF] Teodorescu | Prolonged QRS duration on the resting ECG is associated with sudden death risk in coronary disease, independent of prolonged ventricular repolarization[END_REF][START_REF] Tse | Traditional and novel electrocardiographic conduction and repolarization markers of sudden cardiac death. Eur Eur Pacing Arrhythm Card Electrophysiol J Work Groups Card Pacing Arrhythm Card Cell[END_REF] and are typically assessed by QRS duration and QT. Speckle tracking echocardiography (STE) can identify subtle mechanical alterations of the myocardium affected by fibrosis and electrical conduction disturbances [START_REF] Gorcsan | Echocardiographic assessment of myocardial strain[END_REF]. In healthy subjects, peak longitudinal strain of the different LV segments occurs almost simultaneously with aortic valve closure (AVC), with no or only mild post-systolic shortening and a minimum of mechanical dispersion (MD). In addition, global longitudinal strain (GLS) has become the most robust STE-derived parameter in the prognostic assessment of patients with heart failure and reduced LVEF, dyssynchrony, and SCD [START_REF] Mor-Avi | Current and evolving echocardiographic techniques for the quantitative evaluation of cardiac mechanics: ASE/EAE consensus statement on methodology and indications endorsed by the Japanese Society of Echocardiography[END_REF][START_REF] Kalam | Prognostic implications of global LV dysfunction: a systematic review and meta-analysis of global longitudinal strain and ejection fraction[END_REF]. There is moreover growing evidence supporting the usefulness of temporal deformation parameters (e.g.MD) in predicting outcomes in both ischemic heart disease [START_REF] Hamada | Prediction of Outcomes in Patients with Chronic Ischemic Cardiomyopathy by Layer-Specific Strain Echocardiography: A Proof of Concept[END_REF][START_REF] Huttin | Temporal deformation pattern in acute and late phases of ST-elevation myocardial infarction: incremental value of longitudinal postsystolic strain to assess myocardial viability[END_REF][START_REF] Haugaa | Mechanical Dispersion Assessed by Myocardial Strain in Patients After Myocardial Infarction for Risk Prediction of Ventricular Arrhythmia[END_REF], and heart failure [START_REF] Hasselberg | Left ventricular markers of mortality and ventricular arrhythmias in heart failure patients with cardiac resynchronization therapy[END_REF] as well as in selecting responders to resynchronization therapy [START_REF] Shimamoto | Left Ventricular Mechanical Discoordination in Nonischemic Hearts: Relationship with Left Ventricular Function, Geometry, and Electrical Dyssynchrony[END_REF][START_REF] Sharma | Prognostic implications of left ventricular dyssynchrony for major adverse cardiovascular events in asymptomatic women and men: the Multi-Ethnic Study of Atherosclerosis[END_REF].

While the value of ECG and STE variables has previously been evaluated, their interdependence has not been extensively assessed. Importantly, there is no published threshold in patients without cardiovascular disease and our knowledge of electromechanical coupling is mainly focused on patients with overt cardiac disease. The aims of the present study are to: i) investigate the association between intraventricular conduction delay and delayed repolarization (as assessed by QRS and QTc durations) with that of myocardial deformation profile and dispersion assessed by STE strain analysis, and ii) provide reference values of STE variables according to ECG variables in an initially healthy subjects.

Methods

Study population

The STANISLAS cohort is a single-center familial longitudinal cohort including 1006 families (4295 subjects) from the Nancy region who were recruited in 1993-1995 at the Center for Preventive Medicine [START_REF] Ferreira | Cohort Profile: Rationale and design of the fourth visit of the STANISLAS cohort: a familial longitudinal population-based cohort from the Nancy region of France[END_REF]. The subjects were of French origin and initially healthy, free of acute or chronic diseases (such as hypertension, diabetes, stroke, cancer, etc.). In the latest assessment performed in 2011-2016, a total of 1140 participants (aged 20 to 75 years with no evidence of cardiac disease) who had available echocardiographic images suitable for strain analysis was included in the current analysis (Figure S1). Details of the cohort investigations and loss to follow-up are provided in the Supplemental Digital Content (http://links.lww.com/HJH/A947) and have been published. The research protocol was approved by the local Ethics Committee and all study participants gave written informed consent. Clinical evaluation included information regarding past medical history, symptoms suggestive of any cardiovascular illness, anthropometric measurements (height, weight, and waist circumference), and general physical examination (including blood pressure measurement).

Two-dimensional standard echocardiography

Examinations were performed in the left lateral decubitus position with a commercially available standard ultrasound scanner (Vivid 9, General Electric Medical Systems, Horten, Norway) using a 2.5MHz phased-array transducer (M5S). All echocardiographic and Doppler images were recorded in digital raw-data format and, after performing the centralized anonymization, submitted to the central core laboratory of the Centre d'Investigation Clinique (University of Lorraine, Nancy).

Echocardiographic images were obtained from the 3 standard LV apical views (4-, 2-and 3chambers). All images were acquired at a frame rate of 50 to 70 frames per second for twodimensional views and at least three consecutive cardiac cycles were recorded.

Post-processing analysis of deformation parameters

All strain measurements were performed offline by 2 experienced echocardiographers (OH, SC) with dedicated automated software (Q analysis software, Echo PAC PC version 110.1.0, GE Healthcare). This software enables real-time assessment of strain and strain rate, allowing a very rapid and simple analysis of deformation, by manually placing 3 sample points along the endocardium to define the LV base and the apex on an end-systolic frame and mid-wall myocardial layer strain (equivalent to full-wall) was obtained. Each LV wall was divided into 3 segments (basal, mid and apical), and the quality of tracking was assessed for each myocardial segment. Tracking feasibility in each apical view was rated on visual inspection. All echocardiographic studies with difficult images and/or moderate to poor tracking quality (n=415) were reviewed by a second senior echocardiographer. Segments or views with persistent poor tracking quality were excluded.

Reproducibility analyses showed that inter and intra-observer reproducibility was good, with intraclass correlations >0.70 for all considered parameters [START_REF] Coiro | Reproducibility of echocardiographic assessment of 2D-derived longitudinal strain parameters in a population-based study (the STANISLAS Cohort study)[END_REF].

Strain parameters

LV systolic deformation was assessed using the 3 apical views, with the following parameters derived from strain curves measured for each segment: 1) GLS, corresponding to the maximal absolute value of strain during ejection phase (before AVC), 2) end-systolic strain (ESS), corresponding to the value of strain at end-systole and 3) post-systolic strain (PSS), as the maximal absolute value of strain after AVC during the isovolumic relaxation phase. Post-systolic index (PSI) was automatically calculated by the software as (PSI=100*(PSS-ESS)/PSS), with a value ≥ 8 considered as abnormal.

The time from the beginning of the QRS to the peak early positive (P) systolic (S) and global (G) strain was automatically calculated by the software for each segment from which mean values and standard deviation were subsequently calculated for all patients. A bull's eye representation of these values was obtained (Figure 1). Global values were calculated by averaging all of the values computed at each segmental level from the same frame. To avoid confusion in terminology, a better strain means that the numerical value becomes increasingly negative, while a worst longitudinal strain is observed when LV function deteriorates and GLS becomes less negative. Transmural GLS (GSfullwall) were automatically displayed and calculated from the average of all 17 segments. (Figure 1 & 2)

Electrocardiographic recording

ECGs were recorded at 25 mm/s with an amplification of 0.1 mV/mm in the supine position during quiet respiration and interpreted by a cardiologist. A well-defined protocol for electrode placement was used throughout the study. The RR, PR, and QRS intervals were measured in lead II to the nearest 2 ms from the averaged. Heart rate was calculated using the 1/RR interval. QT interval was measured in lead II or V5 (whichever provided the best delineation of the T wave), with the highest value defined as QTm (QT measured). Using the preceding RR interval, correction of the QT (QTc) interval was obtained by Fridericia's formula (QTcF=QTm/√RR 3 ) [START_REF] Postema | The measurement of the QT interval[END_REF].

Statistical methods

Myocardial deformation parameters were stratified by the quartiles of QRS duration and QTc interval. Mean values of standard echocardiographic parameters, main deformation parameters and temporal deformation parameters were respectively compared between each quartile. Continuous variables are expressed as mean ± SD or median and 25 th -75 th percentiles as appropriate. Histograms and normal-quantile plots were visually inspected to verify the normality of distribution of continuous variables. GLS parameters were normally distributed. Correlation between QRS and QTc delays and STE-derived parameters was evaluated using the Spearman non-parametric method. In order to visualize and to interpret the differences of parameters between Qrs or Qt quartiles separately, radar charts were built using z-scores calculated from mean values of each parameter according to quartiles. P values <0.05 were considered to be statistically significant.

Results

General characteristics

On average, subjects were middle-aged (age 48.6±14.0 years) and overweight (BMI 25.5± 4.2mg/m 2 ). A half (47.7%) were male sex, one fifth (23%) were current smokers and 16.3% had selfreported hypertension or were on anti-hypertensive medication. Mean QRS duration was 91±11ms and mean QTc interval was 408±18ms. Mean values of GSfull-wall was -21.1±2.5%, while mean mechanical dispersion was 33.6±11.7ms.

Association of QRS duration and QT duration with standard echocardiographic parameters

Subjects with longer QRS (increased QRS delay) duration had lower LVEF (Q-QRS=65.9±6.2 vs. Q4-QRS=64.3±6.5; P<0.014), higher LVEDV (Q1-QRS=80.2±20.3 vs. Q4-QRS=104.3±27.7; P<0.001) and higher LV mass index (Q1-QRS=78.4±22.7 vs. Q4-QRS =93.6±29.3; P<0.001). LA volume index was higher in the highest QRS duration quartile (Q1-QRS=21.1±6.4 vs. Q4-QRS=24.0±7.7; P<0.001) (Table 1).

LVEF did not differ across QTc quartiles (P=0.13). LV end-diastolic volume, LV mass index and left atrial (LA) volume index were higher in the highest QTc quartile (all P<0.05). Mean e' was lower in the highest QTc quartile (Q1-QTc=12.1±3.2 vs. Q4-QTc=11.0±3.0, P<0.05) (Table 1).

Association between QRS duration (QRS delay) and myocardial deformation

Strain parameters (GSfull-wall) showed a stepwise decrease across QRS quartiles (all P<0.0001) showing significant worst values in the highest quartile compared to the first two quartiles (P <0.05 for comparison between Q4-Q1 and Q4-Q2) (Table 2); a significant inverse correlation was observed between each of these parameters and QRS delay (P<0.0001) (Table S1). The same significant pattern across quartiles was observed for ESS (Table 2). PSI did not show any significant differences across QRS quartiles. A non-significant pattern across quartiles also emerged with respect to time-to-peak strain and MD parameters (Table 2 & Figure 3).

Association between QT duration (repolarization delay) and myocardial deformation

Contrary to the QRS results, specific strain parameters showed a stepwise increase across QTc quartiles (P<0.0001) with significant better values in the fourth quartile (P <0.05) (Table 2). Significant positive correlations were observed (P<0.0001) (Table S1). Similar results were found for ESS.

All time-to-peak strain parameters significantly increased in a progressive manner across QTc quartiles (all P <0.0001) with significantly better values in the fourth QTc quartile compared to the first two quartiles (P for Q4-Q1 and Q4-Q2 comparisons <0.05) (Table 2); this significant relationship was also confirmed by significant positive correlations (P<0.0001) (Table S1). The same trend was observed for all MD parameters (Table 2

& Table S1)

Normal distribution values by QTc tertiles and sex difference are provided in Table 3.

Discussion

By measuring myocardial deformation patterns in a large population, our results showed that:1) patients with increased QRS duration had lower LVEF and peak systolic strain but no significant relationship with regard to timing and dispersion; 2) patients with longer QTc had a significant delayed time-to-peak strain (about 30ms) associated with an increase in mechanical dispersion. These data thus provide the distribution and range of myocardial deformation parameters according to QTc by sex in a population study.

In addition, in this population-based setting that QRS is not associated with mechanical dispersion, our results suggest that echocardiographic-based dyssynchrony does not largely overlap with ECG-based dyssynchrony. This further supports systematically using imaging to assess dyssynchrony given that the latter is apparently not assessable with routine ECG.

LV mechanics and intraventricular conduction

It is readily acknowledged that cardiac electrical conduction, and cardiomyocyte depolarization culminate in mechanical contraction (excitation-contraction coupling) (1). Mechanical alterations in patients with pathological conduction have been widely described, mainly in heart failure [START_REF] Shimamoto | Left Ventricular Mechanical Discoordination in Nonischemic Hearts: Relationship with Left Ventricular Function, Geometry, and Electrical Dyssynchrony[END_REF][START_REF] Chen | Coupling of ventricular action potential duration and local strain patterns during reverse remodeling in responders and nonresponders to cardiac resynchronization therapy[END_REF], as well as in patients with left ventricular hypertrophy [START_REF] Lin | Electromechanical relationship in hypertrophic cardiomyopathy[END_REF]. An increase in QRS delay leads to mechanical dyssynchrony both between and within the ventricles, a largely-studied phenomenon at the basis of cardiac resynchronization therapy [START_REF] Brignole | ESC Guidelines on cardiac pacing and cardiac resynchronization therapy: the Task Force on cardiac pacing and resynchronization therapy of the European Society of Cardiology (ESC). Developed in collaboration with the European Heart Rhythm Association (EHRA)[END_REF]. In healthy subjects, peak longitudinal strain should occur simultaneously with AVC, with no or only mild post-systolic shortening. Moreover, the peak of strain should also occur simultaneously for the 17 segments with negligible MD. Of particular interest, we did not observe any relationship between QRS duration and MD or with the time-to-peak strain.

Delayed QRS is associated with two pathophysiological phenomena: electrical dyssynchrony and LV hypertrophy [START_REF] Pfeiffer | Biomechanics of cardiac electromechanical coupling and mechanoelectric feedback[END_REF][START_REF] Sharma | Prognostic implications of left ventricular dyssynchrony for major adverse cardiovascular events in asymptomatic women and men: the Multi-Ethnic Study of Atherosclerosis[END_REF][START_REF] Chan | Comparison of the relation between left ventricular anatomy and QRS duration in patients with cardiomyopathy with versus without left bundle branch block[END_REF]. However, in the presence setting, given that our patient cohort did not include patients with QRS over 120ms, it is thus possible to speculate that in patients with normal QRS length, the contribution of QRS to overall mechanical dispersion is minimal.

ECG is the expression of whole heart depolarization and repolarization currents in which their variations in magnitude and timing are reflected by changes in ECG morphology [START_REF] Cardone-Noott | Human ventricular activation sequence and the simulation of the electrocardiographic QRS complex and its variability in healthy and intraventricular block conditions. Eur Eur Pacing Arrhythm Card Electrophysiol J Work Groups Card Pacing Arrhythm Card Cell[END_REF].

We found that patients with increased QRS duration had lower LVEF and peak systolic strain.

As a potential marker of cardiovascular aging and LV/left atrium cavity remodeling, higher LV mass and higher LA volume indices were also observed, which may be related to preclinical diastolic dysfunction [START_REF] Nagueh | Recommendations for the Evaluation of Left Ventricular Diastolic Function by Echocardiography: An Update from the American Society of Echocardiography and the European Association of Cardiovascular Imaging[END_REF]. Interestingly, patients herein with longer QTc showed slightly albeit significantly better peak strain values. This finding, together with the lack of relationship between LVEF and QTc interval, may be explained by the fact that a more prolonged depolarization provides more time for myocardial contraction, which may correlate with better strain values in healthy subjects. The latter contrasts with findings in patients with long QT syndrome in which a small but significantly worst GLS was found when compared with healthy subjects [START_REF] Haugaa | Transmural differences in myocardial contraction in long-QT syndrome: mechanical consequences of ion channel dysfunction[END_REF].

LV mechanics and repolarization delay

According to the present findings, QTc had a more robust relationship with STE-derived temporal parameters than QRS, especially with regard to the delay between electrical activation and peak of deformation. Of note, MD parameters where higher in patients with longer QTc with no relationship with QRS, likely indicating that a longer QRS is likely better correlated with dyssynchrony and not necessarily with prolonged contraction while QTc is the direct measurement of global depolarization duration, including dyssynchrony and prolonged contraction (Online ref .10).

Clinical Implications

Traditionally, isolated LVEF is considered as the predominant parameter for sudden cardiac death (SCD) risk stratification and also as a gatekeeper for resynchronization therapy (2); however, in non-ischemic heart failure, this marker failed to select those who may benefit from implantable cardioverter defibrillator (ICD) in primary prevention of SCD [START_REF] Køber | Defibrillator Implantation in Patients with Nonischemic Systolic Heart Failure[END_REF]. Therefore, novel risk-stratification based on echocardiography is needed (4)(Online ref. 11). In the current study, we found no association between QRS and mechanical dyssynchrony assessed by myocardial deformation imaging (Online ref. [START_REF] Haugaa | Mechanical Dispersion Assessed by Myocardial Strain in Patients After Myocardial Infarction for Risk Prediction of Ventricular Arrhythmia[END_REF], it would appear that the latter has no strong overlap with ECG analysis and may consequently have additional prognostic value on top of QRS data in patients with heart failure. This additional prognostic value would none the less necessitate further evaluation in dedicated studies.

The description of reference values and ranges according to electrical activation of the myocardium in the general population represents a first step in improving awareness of the abnormal LV deformation pattern in myocardial disease. Consequently, we believe the present findings will improve comprehension of electromechanical coupling.

Understanding the interactions between QRS, repolarization and mechanical dyssynchrony will furthermore optimize the application of cardiac resynchronization therapy by improving both patient selection and monitoring [START_REF] Brignole | ESC Guidelines on cardiac pacing and cardiac resynchronization therapy: the Task Force on cardiac pacing and resynchronization therapy of the European Society of Cardiology (ESC). Developed in collaboration with the European Heart Rhythm Association (EHRA)[END_REF]. Indeed, STE provides a practical tool for myocardial mechanics and cardiac synchronism analysis [START_REF] Chen | Coupling of ventricular action potential duration and local strain patterns during reverse remodeling in responders and nonresponders to cardiac resynchronization therapy[END_REF] which is not redundant with routine ECG analysis (as highlighted above).

Study Limitations

Given that QRS durations had little variation between patients (91 ± 10.8ms) in this study, our results cannot be extrapolated to the general population. Subclinical left ventricular hypertrophy may partly explain the results in the higher quartile since impaired systolic and diastolic function was noted in this group.

Since QT interval also includes QRS duration, its length can be overestimated in patients with prolonged QRS (Online ref. [START_REF] Mor-Avi | Current and evolving echocardiographic techniques for the quantitative evaluation of cardiac mechanics: ASE/EAE consensus statement on methodology and indications endorsed by the Japanese Society of Echocardiography[END_REF][START_REF] Hasselberg | Left ventricular markers of mortality and ventricular arrhythmias in heart failure patients with cardiac resynchronization therapy[END_REF][START_REF] Shimamoto | Left Ventricular Mechanical Discoordination in Nonischemic Hearts: Relationship with Left Ventricular Function, Geometry, and Electrical Dyssynchrony[END_REF] and it is therefore important to cautiously interpret the conclusions when a relationship with both QRS and QTc is found in the same direction. In a similar manner, a relationship with QRS and QTc going in opposite directions or only toward one of the latter can be interpreted as a more robust finding.

Conclusion

Physiological variability in QRS and QTc intervals are associated with subtle but significant changes in myocardial deformation patterns. Better longitudinal strain parameters and increased mechanical dyssynchrony relate to longer QTc and increased QRS delay. Integrating the impact of ventricular, repolarization conduction delay, myocardial mechanics and synchrony is likely to improve our understanding of electromechanical coupling.

  Historically, prolonged contraction duration and spatial dispersion of contraction duration were the first alterations identified and initially measured by M mode in patients with long QT syndrome (25,26) (Online ref.1,2). In later studies, tissue Doppler imaging (Online ref.3)(27) and speckle tracking strain analysis (Online ref.4,5)(28), allowed assessing electromechanical changes in this syndrome. In the present study, patients with longer QT had a more prolonged and heterogeneous contraction resulting in higher post-systolic contraction and mechanical dispersion (Online ref.6), thus highlighting the relationship between repolarization duration and dispersion with MD measured by 2D speckle tracking. These findings are in line with a recent study by Sauer et al. showing significant linear association between ECG repolarization heterogeneity (as T wave peak to T wave end) with mechanical dyssynchrony in 82 patients without QT disorders and without significant cardiomyopathy (Online ref.7). The relationship between mechanical dispersion and prolonged QTc may also be associated with myocardial disarrangement as documented by Hurtado-de-Mendoza et al. who showed a correlation between diffuse interstitial fibrosis assessed by T1 relaxation times in cardiac magnetic resonance and QT dispersion in hypertrophic cardiomyopathy (HCM) patients (Online ref.8). QTc prolongation along with MD has moreover been linked to chronically ischemic myocardium and an increased likelihood of coronary artery disease (Online ref.9).

Figure 1 .Figure 2 .

 12 Figure 1. Longitudinal Strain Curves Measured by Speckle Tracking Echocardiography

Figure 3 .

 3 Figure 3. Heat map based on non-parametric correlations of electrical variables, standard

Table 1 . Distribution of Clinical and Standard Echocardiographic Parameters according to QRS and QT Delay Quartiles

 1 § P<0.05 vs. Q2 † P<0.05 vs. Q3. Kruskal Wallis test for nonnormally variable and Bonferonni for nonnormally distributed variables LV: left ventricular; LA: left atrium; TR Vmax: maximum velocity of tricuspid regurgitation.

Table 2 . Distribution of Myocardial Deformation Parameters according to QRS and QT Delay Quartiles GSfull

 2 -wall: global peak longitudinal systolic strain full-wall myocardium; ESS: end-systolic peak longitudinal strain; PSI: postsystolic index (n° of abnormal segments); G: global, S: systolic, P: positive, MD: mechanical dispersion.

Table 3 . Normal Distribution of Myocardial Deformation Parameters according to QTc Tertiles by Sex GSfull

 3 -wall: global peak longitudinal systolic strain full-wall myocardium; ESS: end-systolic peak longitudinal strain; PSI: postsystolic index (n° of abnormal segments); G: global, S: systolic, P: positive, MD: mechanical dispersion

Table 1 : Distribution of Clinical and Standard Echocardiographic Parameters according to QRS and QT Interval Quartiles

 1

	QRS (ms)

* 84.

9 ± 22.9 * 93.6 ± 29.3 * § †

						<0.001	81.2 ± 21.9	84.5 ± 24.5	84.9 ± 22.7	89.5 ± 30.6 *	0.004
	LA volume index, ml/m²	21.1 ± 6.4	22.5 ± 6.9	22.4 ± 6.6	24 ± 7.7 *	<0.001	21.3 ± 6.2	21.4 ± 6.3	22.6 ± 7.2	24.2 ± 7.6 * § †	<0.001
	E/A in mitral valve	1.2 ± 0.42	1.21 ± 0.42	1.18 ± 0.34	1.23 ± 0.41	0.594	1.21 ± 0.4	1.2 ± 0.42	1.19 ± 0.36	1.22 ± 0.43	0.893
	mean e'	11.6 ± 3.1	11.5 ± 3.2	11.7 ± 2.9	11.4 ± 3.1	0.789	12.1 ± 3.2	11.7 ± 3.2	11.4 ± 3	11 ± 3 *	0.001
	TR Vmax, m/s	2.19 ± 0.28	2.14 ± 0.27	2.16 ± 0.31	2.13 ± 0.29	0.281	2.16 ± 0.27	2.15 ± 0.31	2.14 ± 0.29	2.17 ± 0.28	0.807

Table 2 : Distribution of Myocardial Deformation Parameters according to QRS and QT Delay Quartiles

 2

	QRS (ms)

Table 3 . Normal Distribution of Myocardial Deformation Parameters according to QTc Tertiles by Sex

 3

	Male

ACKNOWLEDGEMENTS

Shown in Supplementary material

SOURCES OF FUNDING: