Nina T Odermatt
email: nina.odermatt@mpi-marburg.mpg.de

Moreno Lelli
email: moreno.lelli@unifi.it

Torsten Herrmann
email: torsten.herrmann@ibs.fr

Luciano A Abriata
email: luciano.abriata@epfl.ch

Aleksandre Japaridze
email: a.japaridze@tudelft.nl

Hubert Voilquin
email: hubertvoilquin@gmail.com

Rajkumar Singh
email: rajkumarbiotech@gmail.com

Jérémie Piton
email: jeremie.piton@gmail.com

Lyndon Emsley
email: lyndon.emsley@epfl.ch11giovanni.dietler@epfl.ch

Giovanni Dietler

Stewart T Cole

Prof T Stewart

Cole

Structural and DNA binding properties of mycobacterial integration host factor mIHF

Keywords: Mycobacterial Integration Host Factor, mIHF, Nucleoid Associated Protein, DNA binding, NMR structural determination, atomic force microscopy mIHF, mycobacterial Integration Host Factor, DBS, DNA-binding site, dsDNA, double stranded DNA, NAP, nucleoid associated protein, ssDNA, single stranded DNA

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

Tuberculosis, caused by Mycobacterium tuberculosis, is the leading cause of death by a single infectious agent, claiming over 1.5 million lives each year (WHO, 2016). The pathogen mainly causes lung infections [START_REF] Farer | Extrapulmonary tuberculosis in the United States[END_REF] that are hard to diagnose, difficult to treat and can lead to death when undetected [START_REF] Zumla | Advances in the development of new tuberculosis drugs and treatment regimens[END_REF]. M. tuberculosis can enter a dormant phase [START_REF] Boshoff | Tuberculosis -Metabolism and respiration in the absence of growth[END_REF], characterized by phenotypic drug resistance [START_REF] Gomez | M. tuberculosis persistence, latency, and drug tolerance[END_REF] and has evolved a sophisticated virulence strategy employing type VII secretion systems known as ESX-1 to ESX-5 [START_REF] Gröschel | ESX secretion systems: mycobacterial evolution to counter host immunity[END_REF]. Control of virulence factor expression is crucial and has to be carefully adapted to active growth or the latent phase [START_REF] Zondervan | Regulation of three virulence strategies of Mycobacterium tuberculosis: A success story[END_REF]. Transcription is controlled over interlayered levels involving the DNA supercoiling state and by local transcription factors, such as repressors and activators, sigma factors and nucleoid associated proteins (NAPs) [START_REF] Balleza | Regulation by transcription factors in bacteria: Beyond description[END_REF].

NAPs are small, highly abundant proteins in bacteria with hundreds of binding sites on the chromosome. Their ability to bridge DNA, assisted by supercoiling and macromolecular crowding, leads to compaction of the chromosome [START_REF] Luijsterburg | The architectural role of nucleoid-associated proteins in the organization of bacterial chromatin: a molecular perspective[END_REF] and its organization into 1 Mbp macrodomains [START_REF] Song | Building bridges within the bacterial chromosome[END_REF]. In addition to their architectural role, NAPs act as global gene regulators, often involved in induction of virulence expression and genetic adaptation to environmental cues [START_REF] Dorman | H-NS: a universal regulator for a dynamic genome[END_REF]. In E. coli, at least 12 NAPs have been identified that regulate a vast number of genes [START_REF] Dillon | Bacterial nucleoid-associated proteins, nucleoid structure and gene expression[END_REF]. Although M. tuberculosis is exposed to different stresses during host entry and the accompanying immune response [START_REF] Manganelli | Sigma Factors and Global Gene Regulation in Mycobacterium tuberculosis[END_REF], and can also enter a nonreplicating state, NAPs seem to be underrepresented in this bacterium. To date, only four proteins were confirmed to act as NAPs; Lsr2, HupB, EspR and mIHF. Lsr2 protects M. tuberculosis from reactive oxygen species [START_REF] Colangeli | The multifunctional histone-like protein Lsr2 protects mycobacteria against reactive oxygen intermediates[END_REF], and acts as a dimer [START_REF] Summers | The structure of the oligomerization domain of Lsr2 from Mycobacterium tuberculosis reveals a mechanism for chromosome organization and protection[END_REF] in repressing gene expression [START_REF] Gordon | Lsr2 is a nucleoid-associated protein that targets AT-rich sequences and virulence genes in Mycobacterium tuberculosis[END_REF]. HupB is the only conserved NAP and shows sequence similarity to HU of Gram-negative bacteria. HupB is involved in iron homeostasis, cell wall synthesis and adhesion [START_REF] Pandey | Transcriptional regulation of Mycobacterium tuberculosis hupB gene expression[END_REF]. EspR is essential for the expression of the ESX-1 virulence factors and several other virulence associated genes besides [START_REF] Blasco | Virulence regulator EspR of Mycobacterium tuberculosis is a nucleoid-associated protein[END_REF]. The mycobacterial integration host factor, mIHF, was shown to be essential for bacterial survival and regulates genes involved in lipid metabolism, metabolic pathways, translation and virulence [START_REF] Odermatt | Essential Nucleoid Associated Protein mIHF (Rv1388) Controls Virulence and Housekeeping Genes in Mycobacterium tuberculosis[END_REF].

Most NAPs act as dimers, with each monomer binding to a distinct DNA locus and bringing these loci together upon dimerization, as exemplified by EspR [START_REF] Blasco | Atypical DNA recognition mechanism used by the EspR virulence regulator of Mycobacterium tuberculosis[END_REF], Lsr2 [START_REF] Summers | The structure of the oligomerization domain of Lsr2 from Mycobacterium tuberculosis reveals a mechanism for chromosome organization and protection[END_REF] or HupB [START_REF] Bhowmick | Targeting Mycobacterium tuberculosis nucleoid-associated protein HU with structure-based inhibitors[END_REF] in M. tuberculosis. On the contrary, the mIHF homologue from Streptomyces coelicolor, sIHF, has two separate DNA-binding sites (DBS) and was proposed to bind DNA as a monomer [START_REF] Swiercz | A novel nucleoid-associated protein specific to the actinobacteria[END_REF]. mIHF was shown to bind and bend DNA non-specifically [START_REF] Mishra | Integration Host Factor of Mycobacterium tuberculosis, mIHF, Compacts DNA by a Bending Mechanism[END_REF], but the exact DNA-binding mechanism of mIHF and its structure are unknown. Here, we complement the genetic analysis of mIHF [START_REF] Odermatt | Essential Nucleoid Associated Protein mIHF (Rv1388) Controls Virulence and Housekeeping Genes in Mycobacterium tuberculosis[END_REF] by providing structural insights into its biological function. We determined the solution state structure of mIHF and measured its dynamics by Nuclear Magnetic Resonance (NMR) spectroscopy, and studied its DNA binding properties by NMR, Isothermal Titration Calorimetry (ITC), and Atomic Force Microscopy (AFM).

Methods

Sequence analysis and alignment

Amino acid sequences were obtained from NCBI, aligned with T-Coffee [START_REF] Notredame | T-Coffee: A novel method for fast and accurate multiple sequence alignment[END_REF] and visualized with BioEdit (Thomas A. [START_REF] Hall | BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT[END_REF]. Secondary structure was predicted by PSIPRED [START_REF] Jones | PSIPRED [WWW Document[END_REF].

Cloning, expression and purification of mIHF

Oligonucleotides were synthesized by MICROSYNTH (Switzerland). Primers mihF-exp-F and mihF-exp-R (all primers are listed in Table S1) were used to PCR-amplify the mihF gene beginning at the translational start site. Restriction sites NdeI and BamHI were included in the primers to clone the resulting fragment into pET28a vector to give the expression vector pNO72 producing N-terminally 6x His-tagged mIHF protein. Recombinant His-tagged protein was expressed by growing E. coli BL21 DE3 (ThermoFisher) cells transformed with pNO72 in Luria-Bertani (LB) broth in the presence of 50 µg ml -1 kanamycin. Protein expression was induced at an OD 600 = 0.6 with 500 µM IPTG overnight at 16° C with constant shaking at 180 rpm. Cells were harvested by centrifugation and stored at -80°C until further processing. The pellet was thawed in lysis buffer (50 mM Tris-HCl pH 7.5, 500 mM NaCl, 1 mM Imidazole, 10% (w/v) glycerol and 1% TWEEN 20) supplemented with 1 tablet COMPLETE mini protease inhibitor (Roche) and 10 µl DNase I (Roche) and lysed in a cell disrupter (EmulsiFlex, Avestin, Canada). Insoluble debris was removed after centrifugation and supernatant was loaded on a His-trap column pre-equilibrated with buffer A (500 mM NaCl, 50 mM Tris-HCl pH 7.5).

Increasing concentrations of buffer B (500 mM NaCl, 50 mM Tris-HCl pH 7.5, 500 mM Imidazol) eluted the protein from the column and fractions were visualized by COOMASSIE blue staining of a sodium dodecyl sulphate (SDS) gel INVITROGEN). Pooled fractions containing mIHF were dialyzed into Tris-buffer (150 mM NaCl, 50 mM Tris-HCl pH 7.5) overnight at 4°C. To remove the His-tag, 500 µl thrombin (Sigma) was added to the protein before dialysis. The cut mIHF protein was subsequently purified by passage over Ni-NTA Agarose (Machery-Nagel) and Benzamidine Sepharose (GE Healthcare) to remove uncut protein, His-tags and thrombin. A size-selection purification step was applied to remove unfolded protein and increase sample purity. mIHF was passed over a HiLoad 16/600 column and fractions were pooled and concentrated in centrifugal filters (Ultracel, Amicon) with 3 kDa cut-off. Purified mIHF was stored at -80°C in small aliquots until further use. Concentration of mIHF was measured by Qubit (INVITROGEN).

Primers mihF-R90E-F and mihF-R90E-R were used in the QuickChange II Site-Directed Mutagenesis kit (Agilent) as recommended by the manufacturer. The single amino acid mutation R90E in pNO72 gave rise to pNO84.

For the production of isotope labelled mIHF ( 15 N-and 13 C-ammonium chloride and D-glucose from Cambridge Isotope Laboratories, USA), the exponentially growing culture was resuspended in minimal medium containing 25 mM KH 2 PO 4 , 50 mM Na 2 HPO 4 • 2 H 2 O, 1 g l -1

15 N-labelled ammonium chloride, 10 mM NaCl, 0.2 mM CaCl 2 , 1 mM MgSO 4 , 50 µM iron solution, 10 ml trace metal solution and 3 g l -1 13 C labelled glucose. 1000 x iron solution contained 100 µM citric acid and 50 µM FeCl 3 . 100 x trace metal solution contained 68.2 µM

MnCl 2 • 4 H 2 O, 3.7 µM ZnSO 4 • 7 H 2 O, 0.4 µM CoCl 2 • 6 H 2 O, 0.6 µM CuCl 2 • 2 H 2 O, 1.6 µM H 3 BO 3 , 2.1 µM NiCl 2 • 6 H 2 O, 2.1 µM Na 2 MoO 4 • 2 H 2 O, 1.9 µM Na 2 SeO 3 • 5 H 2 O. E. coli was
grown overnight in minimal medium at 16°C and protein was purified as described above.

Due to the remainder of the His-tag (GSHM), the actual mIHF protein starts with amino acid 5 (underlined). The sequence of the protein used throughout this investigation is thus: GSHMVALPQLTDEQRAAALEKAAAARRARAELKDRLKRGGTNLTQVLKDAESDEVLGKMKVSALLEALP KVGKVKAQEIMTELEIAPTRRLRGLGDRQRKALLEKFGSA.

Circular dichroism (CD) measurement

Circular dichroism of mIHF was measured at a concentration of 10 -20 µM protein in a quartz cuvette of 1 mm path length and analysed using a Jasco J-815 CD spectrometer. Data were converted to mean residual ellipticity and corrected for difference in concentration [START_REF] Greenfield | Using circular dichroism collected as a function of temperature to determine the thermodynamics of protein unfolding and binding interactions[END_REF]. Far was used for protein, and near UV (250 -340 nm) for DNA absorbance. Thermal unfolding was monitored over the whole spectra. Molar ellipticity changes at 222 nm were used to fit a sigmoidal model and calculate the melting temperature (GraphPad Prism). DNA titration was carried out at 20°C with 0, 0.003, 0.006, 0.340, 0.500 and 0.700 equivalents of DNA. A 40 bp dsDNA oligonucleotide (5' -CTGGAGGAGCTGGCAGCAGCGTTTCCGGGTGATGGCTGGT -3') was used. As mIHF does not contain amino acids that absorb at 280 nm, the signal was undisturbed in the transition between near and far UV.

Gel retardation assay

100 ng linear DNA was incubated with increasing amounts of recombinant mIHF (1 pmol -2 nmol) and mIHF R90E (0.1 nmol -2 nmol) for 10 min at room temperature and electrophoresed on a 0.75% agarose gel in 0.5x TBS (25 mM Tris-HCl pH 7.5, 50 mM NaCl).

Post-run staining was performed with GelRed (Biotium) for 20 min at a 1:10´000 ratio of GelRed to 0.5x TBS.

2.5 Nuclear Magnetic Resonance (NMR) spectroscopy and solution structure determination NMR measurements for protein assignment and structure determination were carried out on 1.34 mM uniformly 13 C-and 15 N-labelled mIHF samples in 100 µl, using Shigemi NMR tubes. Protein solution was prepared in 90% H 2 O (50 mM phosphate buffer pH 7.5, 100 mM NaCl) and 10% 2 H 2 O. 0.2% sodium azide was added to prevent sample degradation. The NMR spectra were acquired at 14.1 T, 18.8 T and 23.5 T (600.55, 800.13 and 1000.30 MHz of proton Larmor frequency, respectively) on Bruker Avance III spectrometers, all equipped with a cryogenic cooled probehead. The experiments for determining the 1 H, 13 C and 15 N protein backbone resonance assignment were performed at 14.1 T, using the standard J-UNIO protocol [START_REF] Serrano | The J-UNIO protocol for automated protein structure determination by NMR in solution[END_REF] automated assignment APSY routine. Standard CBCACONH [START_REF] Hiller | APSY-NMR with proteins: Practical aspects and backbone assignment[END_REF], HACACONH [START_REF] Fiorito | Automated resonance assignment of proteins: 6D APSY-NMR[END_REF], and HACANH [START_REF] Hiller | Automated projection spectroscopy (APSY)[END_REF]) APSY experiments were used; the π/2 hard pulses were 11.18 µs 14.45 µs and 37.2 µs for 1 H, 13 C and 15 N, respectively. The number of projections acquired were 56, 56, and 41 with 16, 8, and 8 scans for CBCACONH, HACACONH, and HACANH, respectively, and with a direct acquisition time of 91.7 ms and a recycle delay of 0.9 s. The 1 H-15 N 3D and 1 H-13 C 3D HSQC NOESY spectra used to obtain side-chain resonance assignments and structural restraints for the solution structural determination were performed at 18.8 T and 288 K using standard NMR sequences. The π/2 hard pulses were 10.52 µs 13.6 µs and 35.0 µs for 1 H, 13 C and 15 N, respectively. The mixing times were 100 ms and 120 ms for 1 H-15 N 3D and 1 H-13 C 3D HSQC NOESY, respectively. Up to 1536 x 72 x 100 complex points were acquired for direct 1 H, indirect 15 N and 1 H dimensions, respectively; for 1 H-13 C 3D HSQC NOESY 1088 x 80 x 90 complex points were acquired for direct 1 H, indirect 13 C and 1 H dimensions, respectively. The number of scans were 8, with a direct acquisition time of 84.86 ms and a recycle delay of 1.20 s. Data were processed with 1024 x 128 x 256 complex points matrix and square cosine window functions.

The analysis of the APSY experiments was performed with the software UNIO-MATCH for automated backbone resonance assignment [START_REF] Volk | Automated sequence-specific protein NMR assignment using the memetic algorithm MATCH[END_REF] and UNIO-ATNOS/ASCAN [START_REF] Fiorito | Automated amino acid sidechain NMR assignment of proteins using 13C-and 15N-resolved 3D [1H, 1H]-NOESY[END_REF] for automated side-chain chemical shift assignment. The input for automated NOESY peak picking and NOE assignments with UNIO-ATNOS/CANDID (Herrmann et al., 2002a(Herrmann et al., , 2002b) consisted of the UNIO-MATCH chemical shift assignments for the polypeptide backbone, the UNIO-ATNOS/ASCAN output of side-chain chemical shift assignments, and the NOESY datasets. Table 1 reports NMR structural statistics from NOE assignment and structure calculation using UNIO-ATNOS/CANDID.

DNA titration experiments monitored by NMR

The 15 N-labeled protein was prepared at 200 µM in 50 mM phosphate buffer pH 7.5, 100 mM NaCl, with 10% 2 H 2 O. Titration of this sample with dsDNA (5'-AGCTCGTCAACGCCTT-3', 5 mM stock solution, purchased from MICROSYNTH) was monitored through 1 H, 15 N HSQC spectra collected at 298 K on an AVANCEII-800 MHz spectrometer equipped with a CPTC 1 H, 13 C, 15 N 5 mm cryoprobe. The HSQC spectra were collected with 128 increments in the 15 N dimension. Since broadening was observed for most crosspeaks along the titration, the number of scans was increased throughout the titration as required to keep a workable signal-to-noise ratio (from 8 scans without DNA to 64 scans at 2.2 equivalents of added DNA).

All data for titrations were acquired and processed with a Bruker TOPSPIN 2.0 instrument and analysed with NMRFAM-SPARKY. For data analysis, the backbone 1 H and 15 N assignments derived from the automatic structure calculation procedure were first verified against 3D CBCANH and HN(CA)CO spectra. All assignments except that for Arg27 were confirmed, while assignments for Gly40, Arg92 and Asp96 were also retrieved. After this procedure the only non-proline residues missing H,N assignments are Gly1 through Met4, Leu19, Arg27 through Ala30 and Arg97. The final table of backbone H,N assignments is given in Fig. S1, CBCANH spectra in Fig. S2.

Chemical shift perturbations (CSP) were calculated from weighed differences in 1 H and 15 N chemical shifts (ΔδH and ΔδN, [START_REF] Williamson | Using chemical shift perturbation to characterise ligand binding[END_REF]) as:

𝐶𝑆𝑃 = (∆𝛿𝐻) ! + (∆𝛿𝑁/5) ! 2

15 N-Relaxation data

The backbone dynamics of mIHF was investigated through 15 N-relaxation measurements, acquired on the same sample described above for structure determination at 298 K, corrected for the different temperature. Standard experiments [START_REF] Cavanagh | Protein NMR Spectroscopy, 2nd Edition Principles and Practice[END_REF] were used to measure the longitudinal and transverse 15 N relaxation times ( 1 H-15 N HSQC T 1 and 1 H-15 N HSQC T 2 ) [START_REF] Kay | The measurement of heteronuclear transverse relaxation times in ax3 spin systems via polarization-transfer techniques[END_REF] and the { 1 H}-15 N NOE HSQC at 23.5 T (1.0 GHz of 1 H Larmor frequency). For these experiments, the π/2 hard pulses were 11.93 µs 13.75 µs and 26.5 µs for 1 H, 13 C and 15 N, respectively. The build-up times for the 15 N T 1 were 0.010 s, 0.010 s, 0.100 s, 0.300 s, 0.500 s, 0.700 s, 1.00 s, 1.30 s, 1.60 s, 1.60 s, with two experiments acquired twice. To measure 15 N T 2 a CPMG sequence was used with a τ-π-τ sequence of 450 µs for τ delay and 160 µs for the π-pulse. The CMPG build-up times were 17 ms, 34 ms, 51 ms, 68 ms, 68 ms, 102 ms, 136 ms, 170 ms, 204 ms with one experiment acquired twice. For each buildup point the number of scans were 32 with 1.50 s of recycle delay for the 1 H-15 N HSQC T 1 , and 2.50 s of recycle delay for the 1 H-15 N HSQC T 2 to allow for the probe duty cycle. For the { 1 H}-15 N NOE HSQC 64 scans were acquired for each experiment, with a saturation time of 6.0 s. Each point in the HSQC plane relaxation data was processed with a 2048 x 512 complex points matrix and square cosine windows function for both direct and indirect dimensions.

Relaxation data were analysed with the PROTEIN DYNAMICS CENTER 2.4.6 (2016 Nov/14) software provided by BRUKER BIOSPIN.

Isothermal titration calorimetry

ITC experiments were performed using a Malvern MicroCal PEAQ-ITC calorimeter (Malvern Instruments) at 25°C with 19 injections of substrate addition. ITC cells were loaded with 300 μl of 5-10 μM mIHF in buffer A (10 mM Tris, pH 8.0, 50 mM NaCl). Different length DNA oligonucleotides were used for binding measurements: 7bp ssDNA (3´-CTCACGC-5´, 100-500 μM), 7bp dsDNA (100-250 μM), 16bp ssDNA (5'-AGCTCGTCAACGCCTT-3', 100-250 μM), 16bp dsDNA (100-500 μM), 40bp ssDNA (5´-CTGGAGGAGCTGGCAGCAGCGTTTCCGGGTGATGGCTGGT-3´, 100-250 μM), and 40bp dsDNA (100-500 μM). Data were analysed using Malvern PEAQ-ITC analysis software provided by Malvern instruments.

Atomic force microscopy

The mihF gene and promoter region were obtained by amplification with primers mihFprom-F and mihF-prom-R. Cloning of the product into pTOPO (ThermoFisher) gave rise to pNO75. pNO75 was cut with SpeI and BamHI and the resulting linear fragment was purified with the PCR purification kit (Sigma-Aldrich). I95 cosmid DNA consists of a 42.6 kb region, coordinates 4,292,000 -4,326,000 of the M. tuberculosis chromosome, cloned into pYUB412 and amplified in E. coli [START_REF] Bange | Survival of mice infected with Mycobacterium smegmatis containing large DNA fragments from Mycobacterium tuberculosis[END_REF]. Cosmid DNA was extracted using a Midi Prep kit (Promega) and eluted in TE buffer (10 mM Tris, 1 mM EDTA).

DNA at 0.5 ng µl -1 (final concentration) was mixed with varying concentrations of mIHF, brought to a final volume of 20 µl in buffer II (5 mM MgCl 2 , 5 mM Tris-HCl pH 7.5, 50 mM NaCl, for linear fragment) or in ultrapure water (for cosmid I95). The mixture was deposited on freshly cleaved mica and incubated for 5 min at 25°C. For experiments with I95 cosmid, the mica surface was functionalized with APTES in a separate step prior to DNA deposition as described previously [START_REF] Japaridze | Toward an Effective Control of DNA's Submolecular Conformation on a Surface[END_REF]. AFM images were collected using a MultiMode SPM with a Nanoscope III controller (VEECO Instruments, Santa Barbara, CA, USA) operated in tapping-mode in air. The AFM cantilevers had a spring constant of 5 N/m (Bruker cantilevers, TAP150A) with resonance frequencies ranging between 120 and 160 kHz. All AFM images consist of 512 × 512 pixels with scan frequency ≤1 Hz. Each experiment was performed at least in duplicate and AFM images were obtained at several separate mica locations. Only DNA complexes that were completely visible in an AFM image were considered for statistical analysis. Images were simply flattened using the GWYDDION software (Version 2.22), and no further image processing was carried out [START_REF] Nečas | Gwyddion: an open-source software for SPM data analysis[END_REF].

DNA molecules were traced using DNA TRACE software previously described [START_REF] Mikhaylov | DNA trace: A comprehensive software for polymer image processing[END_REF]. Based on the statistics of tens-to-hundreds of individual molecules the contour length and the radius of gyration were calculated by the software.

LC-MS/MS

Bands of interest were excised from SDS-PAGE gels and in-gel digested with trypsin or chymotrypsin. After gel extraction, samples were dried by vacuum centrifugation and analysed by LC-MS/MS. Peptides were separated on a Dionex Ultimate 3000nano UPLC system coupled in-line to a high-resolution mass spectrometer (ThermoFisher Scientific). A pre-column (ACCLAIM PEPMAP C18 Trap) was used to capture the samples and separation was performed over a 90-min gradient using a capillary column (ACCLAIM PEPMAP C18 RSLC) at 250 nl min -1 . Positive data-dependent acquisition mode was used and precursor peptides with charge > 2 were fragmented. Database search was performed with Proteome Discoverer 1.4

on Mascot, Sequest and MS-Amanda against Tuberculist R27 database ( [START_REF] Lew | TubercuList -10 years after[END_REF].

Met oxidation, Ser-Thr-Tyr phosphorylation and peptide N-acetylation were set as variable modifications while Cys carbamidomethylation was set as fixed modification. Scaffold 4.75 was used for data inspection.

Accession ID

The atomic coordinates and structure of mIHF (Rv1388) have been deposited in the Protein Data Bank under the accession number xxxx.

Results

mIHF DNA binding domain prediction and mutant generation

The 105-amino acid mIHF protein has a calculated molecular weight of 11.5 kDa and contains no tryptophan, tyrosine or cysteine residues, rendering this protein much harder to detect by UV-Vis spectrophotometry. Amino acid sequence alignments (Fig. 1) show that mIHF is highly conserved among Actinobacteria. Conversely, mIHF does not share any sequence similarity with its orthologues from Gram-negative bacteria. Secondary structure prediction revealed mIHF to be an alpha-helical protein, with the first helix (residues 8 -34) predicted to form a coiled coil ( [START_REF] Lupas | Predicting coiled coils from protein sequences[END_REF], Fig. S3). sIHF tuberculosis [START_REF] Swiercz | A novel nucleoid-associated protein specific to the actinobacteria[END_REF].

Alignment of mIHF with sIHF showed high homology levels between the two sequences, with 62% identity and 79% similarity ( [START_REF] Altschul | Basic local alignment search tool[END_REF], BLOSUM62 matrix). These data suggested that the two proteins might have a very similar structure. We therefore exploited the sIHF structure with the predicted DBSs [START_REF] Swiercz | A novel nucleoid-associated protein specific to the actinobacteria[END_REF] to identify the putative DBSs of mIHF. Overall, both domains were conserved but several interesting differences between sIHF and mIHF were noticed. For instance, positions G67 and G80 in sIHF were occupied by the positively charged lysine and the negatively charged glutamic acid residues, respectively, in mIHF. These changes might reflect the adaptation to specific DNA binding patterns in the respective genomes or modify sequence specificity. Furthermore, several positively charged arginine residues are present in the mIHF and sIHF DBS and are presumably responsible for the direct contact with the negatively charged DNA backbone. We chose R90, R92 and G93 conserved in both proteins, as targets of site-directed mutagenesis to investigate the DNA binding properties of the mutant proteins. Arginine was mutated to glutamic acid (R90E, R92E), and glycine to tryptophan (G93W). The mutant proteins containing R92E or G93W

were not stable after purification, while mIHF R90E was as stable as mIHF wild type (wt) and therefore used for future experiments.

mIHF is a highly soluble, alpha helical DNA-binding protein

The mIHF protein was purified to homogeneity with a concentration of 60 mg ml -1 and displayed elevated stability at room temperature. The structural features of mIHF were characterized via circular dichroism (CD) spectroscopy. Fig. 2a illustrates the mean residue ellipticity of mIHF wt and its R90E mutant. Both proteins showed a folded structure dominated by alpha helices, as reflected in the two pronounced negative peaks at 208 and 222 nm. mIHF R90E displayed a minor deviation from mHIF wt, suggesting that their overall secondary structures are very similar but not completely identical. The CD signal at 222 nm, as a function of temperature, was fitted to a

a) b) c) d) e) f)
sigmoidal non-linear model and the melting temperature was calculated as 52°C (Fig. 2b).

The protein did not refold into its alpha helical structure at 20°C after heating to 95°C. The melting temperature increased to 67°C in the presence of DNA, suggesting higher stability of mIHF when bound to DNA (Fig. 2b). To screen for conformational changes in the mIHF protein upon DNA binding, the wt and mutant proteins were incubated with very low (0.003 and 0.006 equivalents) and increasingly higher amounts (0.34, 0.5, 0.75 equivalents) of a 40 bp dsDNA oligonucleotide. The minima at 208 nm and at 222 nm became less pronounced in mIHF wt at low amounts of DNA, indicating a small conformational change with no emergence of a different secondary structure. At higher DNA/protein ratios, a deeper minima at 208 nm appeared (Fig. 2c and zoom in Fig. 2d). mIHF R90E showed a similar pattern with a decrease in both 208 nm and 220 nm minima, but no subsequent increase at high DNA concentrations. To assess how the mIHF protein shifted linear DNA, a 700 bp DNA fragment was incubated with mIHF wt and mIHF R90E in a gel retardation assay. While mIHF wt shifted the DNA effectively, mIHF R90E had a markedly lower retardation effect (Fig. 3). This confirmed that R90 is important for DNA binding by mIHF. The minor DNA shift caused by mIHF R90E suggested that the mutant protein is still capable of binding DNA, confirming the CD results, but with a decreased DNA binding strength compared to mIHF wt. Thus, mIHF R90E represented a suitable control for other DNA-mIHF interaction experiments.

Solution structure of mIHF and protein dynamics

In order to define the exact binding mechanism and interaction with DNA, structure determination of mIHF was carried out in solution by NMR spectroscopy. High-dimensional automated projection spectroscopy (APSY) spectra for backbone assignment and 3D 15 N/ 13 C-HSQC-NOESY spectroscopy to obtain side-chain resonance assignments and conformational restraints were acquired (Table 1). The backbone and total atom assignments were completed to 90.24% and 76.71%,

respectively. An additional four amino acids (GSHM) remained after removal of the His-tag;

the mIHF protein therefore starts at residue number 5. Note that the NMR statistics (Table 1) concern the entire construct from residues 1 to 109, and the residue numbering in the following section refers to the protein construct used for the NMR study. mIHF was identified as a globular protein (Fig. 4a) for residues 12 -107 with an RMSD (root mean square deviation of atomic positions) of 1.45 Å or 1.11 Å for residues 16 -107, excluding the first four N-terminal helical residues which show increased conformational spread. The Nterminal domain of the mIHF protein (residues V5-T11) did not form a defined secondary structure and was disordered (Fig. 4b) which is consistent with the NMR chemical shift-based secondary prediction. Overall, the NMR structure of mIHF resembled closely the sIHF X-ray structure with an overlay backbone RMSD of 1.62 Å and a similar topology with six α-helices. The flexible Nterminal end of mIHF is constituted by residues V5 to R15 followed by the first α-helix (α1) comprising A16 to R38. This helix was predicted to form a coiled coil and might be the site of dimerization for mIHF. A minor kink at residues A28 and R29 caused a slight divergence from an optimal helix. α1 is connected by a five-residue loop to α2, a short helix from L47 to E54.

α3 only spans one turn from K57 to V61 and is followed by α4 from L64 -K70, α5 (Q77 -P87) and finally α6 (A101 -A109) (Fig. 4b).

Backbone assignments by UNIO-MATCH were completed to 90.24%. To probe the internal dynamics of mIHF and its aggregation state, three relaxation parameters for each backbone amide were characterized: the longitudinal and transverse 15 N relaxation rates R 1 =1/T 1 and R 2 =1/T 2 , respectively, as well as the steady state heteronuclear Overhauser effect ({ 1 H}-15 N NOE). 15 N relaxation data are related to the correlation time of the reorientation motion of the H-N bonds due to protein dynamic in solution. For globular proteins in solution this is primarily due to the molecular rotational correlation time, which in turn is related to the molecular hydrodynamic radius. We can derive the molecular rotational correlation time from the T 1 /T 2 ratio for the 44 most rigid residues (with { 1 H}-15 N NOE >0.65, and neglecting the internal dynamics) that provided a rotational correlation time of 6.9 ± 1.7 ns (at 298 K), which is in line with a molecular weight of 10-12 kDa. In the supporting information the measured values of 15 N R 1 , R 2 and { 1 H}-15 N NOE were reported in table S2 and plotted in Fig.

S4.

3.4 mIHF binds to double-stranded DNA and has two distinct DNA-binding domains

The ligand binding affinity of purified mIHF protein with dsDNA was determined by isothermal titration calorimetry (ITC). mIHF can bind dsDNA of different lengths, with affinity constants in the low micromolar range. The dissociation constant (K D ) values for 7, 16 and 40 bp dsDNA were 6.08, 2.39 and 9.16 μM, respectively (Fig. S5 a, b, c). Binding was only observed for dsDNA, but not for ssDNA oligonucleotides of the same lengths (Fig. S5 d,e,f).

To investigate the interaction of mIHF with dsDNA in structural terms, we performed a titration of 200 µM 15 N-labelled protein with 16 bp dsDNA and monitored the effects of DNA addition on mIHF by solution NMR (Figure 5). Specifically, we titrated a dsDNA molecule of sequence 5'-AGCTCGTCAACGCCTT-3' into 15 N-labelled protein and acquired 1 H-15 N HSQC spectra at 0, 0.25, 0.50, 1.0 and 2.20 equivalents (dsDNA/protein ratios). Analysis of the titration data relied on the backbone 1 H and 15 N assignments obtained from the automated assignment and structure calculation procedures, which we verified and completed by inspecting CBCANH and HN(CA)CO 3D spectra (See fig. S1 / S2 for confirmation of assignments). Of the non-proline residues, only the four N-terminal residues plus L19, R27 through A30 and R97 lack backbone H and N shifts in our final assignment, and hence were not included in the analysis.

Effects on mIHF HSQC spectra were evident already at 0.25 equivalents, as we observed the shifting and broadening of several crosspeaks. Figure 5A illustrates this on a region of the HSQC spectrum, with well-resolved crosspeaks, that exemplifies K70

broadening beyond detection and several other crosspeaks experiencing large shifts.

Resonance broadening suggests exchange kinetics in the millisecond-microsecond timescale, which is consistent with an affinity in the micromolar range as determined by ITC. To better define the chemical shift perturbations, we continued the titration recording spectra at 0.50, 1.0 and 2.2 equivalents, increasing the number of scans as required to compensate for the broader signals. Figure 5B presents a plot of intensity ratios at 0.25 equivalent (where effects are already strong but most signals are still intense enough for quantification) and Figure 5C presents the chemical shift perturbations at 2.2 equivalents (where they are largest). Figure 5D in turn maps these intensity drops and chemical shift perturbations onto the mIHF structure. We observe clustering of the residues most affected by intensity loss or shifts in two sites, which is consistent with dsDNA bound to the two sites identified in the X-ray structure of sIHF bound to DNA [START_REF] Swiercz | A novel nucleoid-associated protein specific to the actinobacteria[END_REF]. This is relevant because it means that the two sites identified for sIHF across the crystallographic unit cell are not merely an artefact of crystallization, but actual binding sites functional in solution. Residues K33, K37, E66, L68, K70, G72, K75 and A76 are all impacted by broadening (red in Figure 5D) and thus correspond to DBS-I of sIHF, together with K73, which experiences a large shift (blue).

Residues K60, S62, L68, V71, M80, T81, T88, R89, L91, R92, G93, D96, Q98, were largely affected by chemical shift perturbation (blue in Figure 5D) and correspond to DBS-II proposed for sIHF together with R90 which experiences broadening (red). Qualitatively, it appears that binding to DBS-II proceeds at a faster exchange timescale, i.e. is of lower affinity than DBS-I.

Finally, from the results of our titration, we also observed some residues undergoing chemical shift perturbations and broadening outside of the DNA interaction sites, especially on one face of the N-terminal helix: A17, K21, A23, A24. Notably, this helix has some predicted coiled coil propensity (fig. S3), which could imply weak dimerization upon DNA binding.

mIHF introduces defined loops into linear DNA fragments

To study mIHF binding to DNA, purified mIHF protein was incubated with DNA substrates of various length, topology and GC content, and imaged using AFM. The first substrate tested was a 1 kb linear fragment containing the mihF gene and the upstream promoter region with a GC content of 65%, representing the mean GC content in M. tuberculosis. The DNA fragment alone showed a looping probability of 10%, which increased upon addition of mIHF and reached 53% when the maximum amount was added (Table S3, Fig. 6). The size of the loops was well defined with a mean length of 110 ± 10 nm, which corresponds to ca. 300 bp (Table S4). When the same experiment was repeated with mIHF R90E, loops in the linear DNA were observed as well but with a lower ratio compared to the wt protein. Only 28% looped fragments were detected upon binding of mIHF R90E at the same maximal concentration (Table S3, Fig. 6b). The average loop length was similar for mIHF and mIHF R90E with 110.5 ± 10 nm and 106.75 ± 10 nm, respectively (Table S4). These results indicated that mIHF R90E did exhibit DNA binding activity, but to a lower extent compared to the wt protein, confirming the gel retardation assay results. Naked DNA had approximately 50% left-and right-handed loops, corresponding to positive and negative supercoiling, respectively. mIHF wt, as well as mIHF R90E, increased the number of lefthanded loops in a concentration-dependent manner to 76% for the wt protein and to 68%

for the mutant (Fig. 6c, Table S5).

Large cosmids unwind and compact in an mIHF concentration dependent manner

The M. tuberculosis chromosome is a circular, supercoiled macromolecule of 4.4 Mb in length. Topology and size are not correctly represented by short linear DNA fragments or plasmids. Furthermore, bacterial chromatin is organized into microdomains of ca. 10 kb [START_REF] Dame | Bacterial chromatin: converging views at different scales[END_REF], a compartmentalization that cannot be mimicked by plasmids.

Therefore, we used the 42.6 kb I95 supercoiled cosmid [START_REF] Bange | Survival of mice infected with Mycobacterium smegmatis containing large DNA fragments from Mycobacterium tuberculosis[END_REF] to reproduce more closely the features of the genomic DNA. Supercoiled I95 forms highly ordered structures, the so-called hyperplectonemes, only present in DNA molecules above 30 kb in size [START_REF] Japaridze | Hyperplectonemes: A Higher Order Compact and Dynamic DNA Self-Organization[END_REF]. Incubation of I95 with mIHF wt at low concentrations (0.6 mIHF bp -1 ) led to unwinding of the cosmid from the hyperplectonemic form and reduced its complexity (Fig. 7). At higher mIHF ratios (1.2 to 3.6 mIHF bp -1 ), the cosmid collapsed and condensed again, as illustrated by the decrease in radius of gyration (smallest circle that contains the DNA polymer, Fig. 7b, Table S6). Re-compaction of the cosmid at high mIHF ratios only reached first-order plectonemes, but did not lead to the formation of hyperplectonemes like those observed in the absence of protein. Similar to the loops observed on the linear DNA fragment, mIHF also introduced mostly left-handed loops in the cosmid, as illustrated in Fig. 7c andd. The R90E mutant simplified the topology of I95, but no collapse was visible at higher ratios of the protein, in agreement with the reduced DNA-binding observed with linear DNA (Fig. S6).

Discussion and conclusion

The importance of mIHF in M. tuberculosis physiology was suggested by Schubert and colleagues, whose mass spectrometry experiments predicted the protein to be among the ten most abundant proteins in the bacterium [START_REF] Schubert | Absolute Proteome Composition and Dynamics during Dormancy and Resuscitation of Mycobacterium tuberculosis[END_REF]. This initial indication was confirmed by our genetic studies, where we proved the essentiality of the gene and the role of the protein in controlling transcription of a variety of virulence factors and housekeeping genes [START_REF] Odermatt | Essential Nucleoid Associated Protein mIHF (Rv1388) Controls Virulence and Housekeeping Genes in Mycobacterium tuberculosis[END_REF]. Here, we investigated the structure of mIHF and its DNA binding properties.

High similarity to the homologue sIHF from S. coelicolor was already evident from the primary structure alignment and further validated by the NMR structure, which differed for only 1.3 Å compared to the sIHF crystal structure. mIHF formed a globular protein with a protruding α-helix at the N-terminal end and the core of the protein consisted of five short α-helices.

ITC confirmed a K D in the low micromolar range, corresponding to a medium-high affinity of the transcription factor for dsDNA. Further, the difference in K D for different substrate lengths ranged only four-fold, suggesting no preference for DNA length. No binding was observed to ssDNA, and similar results were obtained by Swiercz et al., who showed a much higher affinity of sIHF for dsDNA than for ssDNA [START_REF] Swiercz | A novel nucleoid-associated protein specific to the actinobacteria[END_REF].

ITC and NMR titrations of mIHF with DNA, and temperature stabilization of mIHF in presence of DNA, all clearly indicate that mIHF binds specifically double stranded DNA molecules. In particular, mapping of binding sites by NMR spectroscopy indicates binding on two sites that are consistent with those identified in the X-ray structure of sIHF. High-resolution AFM images showed that incubation of mIHF with a 1 kb linear DNA fragment introduced loops of a defined size. Such DNA loops are commonly associated with regulation of gene expression in bacteria, best described in the lac operon of E. coli, where LacI activates transcription by looping the promoter region and by constraining supercoils [START_REF] Fulcrand | DNA supercoiling, a critical signal regulating the basal expression of the lac operon in Escherichia coli[END_REF]. Although mIHF does appear to bind DNA in a sequence-independent manner, it affects the expression of many genes in M. tuberculosis [START_REF] Odermatt | Essential Nucleoid Associated Protein mIHF (Rv1388) Controls Virulence and Housekeeping Genes in Mycobacterium tuberculosis[END_REF]. The mIHF protein might therefore have a direct effect on gene expression by looping the DNA and by bringing other regulatory elements closer to the target gene. Support for this is provided by the large number of loci occupied by mIHF and other NAPs, such as EspR [START_REF] Odermatt | Essential Nucleoid Associated Protein mIHF (Rv1388) Controls Virulence and Housekeeping Genes in Mycobacterium tuberculosis[END_REF]. 76% of the loops observed at the highest protein:DNA ratio (1.22 mIHF bp -1 ) were of left-handed orientation.

Naturally, DNA in bacterial cells is tightly packed in negative supercoils as exemplified by the I95 cosmid. Introduction of left-handed loops by mIHF relaxed and opened the cosmid, thus making it potentially accessible to replication and transcription machineries. Indeed, mIHF was found to act mainly as a transcriptional activator [START_REF] Odermatt | Essential Nucleoid Associated Protein mIHF (Rv1388) Controls Virulence and Housekeeping Genes in Mycobacterium tuberculosis[END_REF]. This pattern is similar to that caused by E. coli HU on cosmid I95, where ordered, regular, left-handed loops were observed [START_REF] Japaridze | Hyperplectonemes: A Higher Order Compact and Dynamic DNA Self-Organization[END_REF]. Similarly to mIHF, E. coli HU binds without sequence specificity to DNA, but prefers gapped and cruciform DNA [START_REF] Pinson | Differential binding of the Escherichia coli HU, homodimeric forms and heterodimeric form to linear, gapped and cruciform DNA[END_REF]. On the contrary, M. tuberculosis EspR introduced loops and additional twists into the DNA [START_REF] Blasco | Virulence regulator EspR of Mycobacterium tuberculosis is a nucleoid-associated protein[END_REF]. Therefore, it appears that mIHF and EspR have divergent effects on the DNA topology of the tubercle bacillus.

In conclusion, we showed here that mIHF has a globular structure, comprising six alpha helices, that is similar overall to that of sIHF despite the divergent sequences. The protein has two DNA binding sites that are functional in solution. DNA binding plays an important role in shaping the DNA structure by introducing left-handed loops. We propose these DNA binding capabilities could possibly be involved in promoting transcription.

Acknowledgments

We would like to thank the mass spectrometry facility at EPFL for their technical support and Stefanie Boy-Röttger for technical support. This work was supported by the Swiss National Science Foundation (grant 31003A-162641 to STC) and by H2020-INFRAIA iNEXT (grant 653706) for supporting the NMR instrument time in the Lyon CRMN platform.

Author contributions

Fig. 1 .

 1 Fig. 1. Amino acid sequence alignment of IHF proteins. Amino acids are similarity coded using the BLOSUM62 matrix at 80% cutoff. mIHF from M. tuberculosis (M.tb) is aligned with protein sequences from M. marinum (M.ma), M. ulcerans (M.u), M. haemophilum (M.h), M. leprae (M.l) and sIHF from S. coelicolor (sIHF, SCO1480).

  contains a helix-two turns-helix (H2TH) motif, where the two turns between the helices contact DNA at DBS-I of the protein. The second DBS lies on the opposite side of the protein, in a turn connecting two helices[START_REF] Swiercz | A novel nucleoid-associated protein specific to the actinobacteria[END_REF] and the proposed DNA binding mode suggests a new category of NAPs, characterized by two distinct DNA binding domains in the same protein rather than by dimerization of two proteins with one DBS each. sIHF is thus unique and not related to other NAPs like IHF or H-NS from E. coli or Lsr2 from M.

Fig. 2 .

 2 Fig. 2. Circular dichroism (CD) spectra of mIHF protein mutants. a) CD spectra of mIHF wt and mIHF R90E. b) Thermal unfolding of mIHF wt in absence and presence of dsDNA at 222 nm. Dotted lines indicate melting temperature with (red) and without (black) DNA. c) CD spectra of mIHF wt with DNA titration and zoom of minima representing alpha helices (d). Different DNA equivalents were used as indicated. e) CD spectra of mIHF R90E with DNA titration and zoom (f), as for mIHF wt.

Fig. 3 .

 3 Fig. 3. Gel retardation assay obtained upon incubation of a 700 bp DNA fragment with mIHF or with mIHF R90E. DNA was incubated with increasing concentrations of mIHF as detailed in the image and run on a 0.75% agarose gel. Post-run staining was performed with GelRed. Ladder: 100 bp molecular weight marker.

Fig. 4

 4 Fig. 4 Structure of mIHF. a) Ribbon drawing of mIHF (90° rotated around z-axis). b) Overlay of the 20 lowest energy structures of mIHF. c) Surface charge of mIHF, blue representing positively and red negatively charged surface of mIHF. DNA (red) is superimposed from PDPB ID 4ITQ after alignment of the sIHF protein to mIHF.

Fig. 5 .

 5 Fig. 5. Titration of double-strand DNA onto mIHF, monitored by 1 H, 15 N HSQC spectra. a) Titration of 15 N-mIHF with DNA, followed through NMR spectroscopy at 800 MHz 1 H frequency. Overlay of 1 H-15 N HSQC spectra without DNA (red) and with 0.25 (orange), 0.50 (green), 1 (purple) and 2.2 equivalents of 16 bp dsDNA (blue). b) Residue-wise ratio of crosspeak intensities at 0.25 equivalents (DNA/protein) relative to the intensities in the spectrum of the free protein. (c) Residue-wise 1 H, 15 N chemical shift perturbations at 2.2 DNA/protein equivalents relative to the free protein. In panels b and c, residues labeled with an asterisk lack H,N assignments or are prolines. (d) Mapping of residues displaying a sharp drop in the intensity ratios at 0.25 equivalents (red spheres, I/I0 < 0.6) and of residues experiencing large chemical shift perturbations (blue spheres, CSP > 0.08 ppm) on a cartoon representation of mIHF NMR structure. Green spheres are for residues that experienced both large intensity drop and large chemical shift perturbations. The two dsDNA molecules are overlaid from PDB ID 4ITQ after alignment of the sIHF protein to mIHF. Cartoon portions in black correspond to residues missing H,N assignments, marked with an asterisk in panels (b) and (c). A 3D interactive version of panel b is available at http://lucianoabriata.altervista.org/papersdata/mihfdnatitration.html

Fig. 6 .

 6 Fig. 6. Atomic force microscopy experiments. A 1 kb linear DNA fragment harbouring the mihF gene and its promoter region was incubated with the mIHF protein. a) Image of DNA fragments and zoomed images of looped structures at 0.21 mIHF bp -1 . Scale bar represents 500 nm. b) and c) Graphs representing the mean percentage ± standard deviation of looped fragments observed upon addition of increasing concentrations of mIHF wt and mIHF R90E and the percentage of left-handed loops.

Fig. 7 .

 7 Fig. 7. Atomic force microscopy experiments with cosmid I95. a) Cosmid I95 alone and incubated with increasing concentrations of mIHF wt. The lower panels are zoomed images of the squares defined in the upper panels. Scale bar represents 1 μm and mIHF bp -1 ratios are indicated above images. b) Radius of gyration measured for cosmids in the presence of mIHF. Bars represent mean ± standard deviation. c) Cosmid I95 incubated with 1.2 mIHF bp -1 , zoom to an example structure and handedness of loops ("-" = left-handed, "+" = right-handed, (d)).

  NTO conceived the experiments, carried out the wet lab experiments, ML and TH determined the NMR structure of mIHF, LAA verified resonance assignments and performed and analysed the NMR titration, AJ and HV did the AFM studies, RS performed the ICT experiment, JP suggested the mutations, LE, GV and STC supervised the study. NTO, LAA and STC wrote the manuscript.

Table 1 . NMR structural statistics. NMR constraints

 1

	Distance


A Pineau

G Guillemot

G Reinhart

G Regula

N Mangelinck-Noël

Ch-A Gandin

3D cellular automaton modelling of silicon crystallization including grains in twin relationship

Production of silicon for solar cells in photovoltaic systems is mainly based on directional casting processes. Twin nucleation is favoured during silicon growth due to the lowlevel twin energy of formation. As a consequence, in all solidification process, a large amount of grain boundaries (GB) are in twin relationship. A 3D cellular automaton (CA) model has been recently developed for the growth of multi-crystalline silicon including facet formation and nucleation of new grains in twin relationship. Activation of facets is based on an undercooling parameter assigned to each grain and for each of the <111> crystal directions. The model also considers nucleation and growth of grains on <111> facets corresponding to Σ3 twin relationships between twin grains. This model is first applied to comparison with experimental observations. It is found that impingement of growing grains that nucleated in Σ3 twin relationships meet during growth and form Σ3, Σ9 and Σ27 GB, in good agreement with experimental observations. Finally, the model is applied at a larger scale to generate grain structures representative of industrial practice. While quantitative experimental data is missing for comparison at such scale, the model is promising and its implementation in heat and mass transfer models should be considered for assistance to production of silicon for solar cells dedicated to photovoltaic systems.

Introduction

Industrial photovoltaic panels are mainly based on the use of silicon material. Crystalline silicon has the largest market share due to: the availability of the raw material, non-toxicity, efficiency that could reach 15% and a high level of technology readiness [START_REF] Lacerda | Diversity in solar photovoltaic energy: implications for innovation and policy Renew[END_REF]. The photoelectric properties of the solar cells are strongly linked to the production chain, including purity of the raw material, crystallization, wafering, cell processing, and module integration. During crystallization, the grain structure and associated defects are generated, defining density, crystallographic orientations and twin-relationships. Similarly, impurities are developed and dislocations are formed in material. The crystallization features need to be controlled in order to optimize final electrical properties because they can act as recombination centres of charge carriers, thus lowering the efficiency of the photoelectric panels [START_REF] Trempa | Defect formation induced by seed-joints during directional solidification of quasi-mono-crystalline silicon ingots[END_REF][START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 〈110〉 seeded growth[END_REF][START_REF] Kakimoto | Crystallization of Silicon by a Directional Solidification Method, Crystal Growth of Si for Solar Cells[END_REF][START_REF] Ouaddah | Role of impurities in silicon solidification and electrical properties studied by complementary in situ and ex situ methods[END_REF].

During the manufacturing of silicon ingots in crucibles, twin grain boundaries are commonly observed in any solidification process and strongly influence the final grain structure. Voigt et al. [START_REF] Voigt | Proc 14th Europ. PV Solar Energy Conf[END_REF][START_REF] Voigt | Grain orientation and grain boundaries in cast multicrystalline silicon[END_REF] analysed the orientation of each grain on a 5×5 cm 2 wafer. They observed that most of the grains are in twin relationship and more than 50% of the grain boundaries correspond to Σ3, Σ9, Σ27 and Σ81 twins boundaries. The most frequent twins in silicon are of Σ3 type and can be described as an odd number of 60° rotation of the crystal lattice around a <111> direction. New grains in Σ3 relationship nucleate on {111} plane [START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 〈110〉 seeded growth[END_REF], their formation being enhanced by the low energy needed for twinning process [START_REF] Bristowe | Properties of Crystalline Si[END_REF]. According to Jackson's criterion [START_REF] Jackson | Constitutional supercooling surface roughening[END_REF], facets may correspond to the single {111} crystallographic planes in silicon material. Being coherent boundaries, Σ3 twins exhibit no dangling bonds and do not directly lower the final efficiency of solar cells [START_REF] Duffar | Comprehensive review on grain and twin structures in bulk photovoltaic silicon[END_REF]. However, the nucleation of a grain in Σ3 twin relationship can be at the origin of the formation of higher order incoherent, symmetric or asymmetric twin boundaries such as Σ9 or Σ27 [START_REF] Ribéri-Béridot | On the impact of twinning on the formation of the grain structure of multicrystalline silicon for photovoltaic applications during directional solidification[END_REF] and hence responsible for part of the dislocation emission during growth [START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 〈110〉 seeded growth[END_REF][START_REF] Ryningen | Growth of dislocation clusters during directional solidification of multicrystalline silicon ingot[END_REF] and for impurity segregation thus degrading the photoelectric properties.

CA aims at modelling complex phenomena through the use of local simple laws. This approach was firstly proposed by Gandin and Rappaz [START_REF] Ch | A coupled finite element-cellular automaton model for the prediction of dendritic grain structures in solidification processses[END_REF] to model solidification processes. The development of grain envelopes is simulated rather than the full dendritic microstructure, with the advantage to be applicable to large simulation domains. The present work aims at modelling silicon growth upon directional solidification using a 3D CA method. The proposed model is firstly introduced. It is based on phenomenological and geometrical considerations, so as to retrieve {111} facet growth and nucleation of grains in twin relationship as commonly observed experimentally. Finally, the model is applied at a macroscopic scale to highlight the influence of facet growth and new twinned grains nucleation on the final ingot grain structure.

Modelling

The description of the CA algorithm is presented in [START_REF] Pineau | Cellular automaton modelling of silicon crystallization with grains in twin relationship[END_REF] with details and associated schematics. The strategies implemented in the CA model to allow modelling of polycrystalline silicon growth are here briefly summarized. The model used in the present work is divided into three steps: (i) {111} facet activation and growth, (ii) grain in twin relationship nucleation and (iii) twin growth. In the following subsections, index μ refers to a liquid cell.

Facet activation and growth

The standard growth algorithm developed for cubic materials is used [START_REF] Pineau | Cellular automaton modelling of silicon crystallization with grains in twin relationship[END_REF]. However, dedicated developments are introduced to model faceted growth. ).

This stochastic approach is necessary to reproduce the size distribution between twin grain boundaries.

The mean value and standard deviation parameters are adjusted based on comparison with experimental observations.

Facet growth.

Facets of a grain 𝑔 being described by its Cartesian equations 𝐷 𝑔 {111} : 𝑎 𝑔 {111} 𝑥 + 𝑏 𝑔 {111} 𝑦 + 𝑐 𝑔 {111} 𝑧 + 𝑑 𝑔 {111} = 0, their growth is simply modeled by the time integration of the crystal growth velocity of the {111} facets 𝒗 𝑐 {111} defined along the normal 𝒏 {111} . In practice it is sufficient to update 𝑑 𝑔 {111} with time using

𝑑 𝑔 {111} = 𝑣 𝑐 {111} (𝑡 -𝑡 0 ) + 𝑑 𝑔,𝑡 0 {111} (1)
where 𝑣 𝑐 {111} is taken constant following measurements available elsewhere [START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF]. The CA standard algorithm is modified [START_REF] Pineau | Cellular automaton modelling of silicon crystallization with grains in twin relationship[END_REF] and referred as conditioned capture algorithm. This is achieved by conditioning the capture of a cell  by a neighboring growing cell belonging to grain 𝑔 considering its location with respect to the facets described by their Cartesian equations. [START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 〈110〉 seeded growth[END_REF][START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF]. Upon capture of a liquid cell  by its neighboring growing cell , location of cell  is tested to know if it is located on a sample edge or as part of the bulk. For an edge cell (respectively a bulk cell), if the test Δ𝑇 𝜇 > Δ𝑇 𝑁,𝑇 𝐸 (respectively Δ𝑇 𝜇 > Δ𝑇 𝑁,𝑇 𝐵 ), is verified, a nucleation event takes place. The Euler angles of the new grain formed in cell 𝜇 correspond to the crystallographic orientations of the capturing grain to which a Σ3 twin rotation (60°) around the normal to the capturing facet is applied.

Twins growth

Twin boundaries between grains are locked during growth [START_REF] Tsoutsouva | In situ investigation of the structural defect generation and evolution during the directional solidification of 〈110〉 seeded growth[END_REF][START_REF] Stamelou | {111} facet growth laws and grain competition during silicon crystallization[END_REF]. To maintain the twin boundary while pursuing growth of the solid-liquid interfaces, the conditioned capture algorithm is used. The twin grain boundary remains static in space and develops as a result of growth of the grains on each side of the corresponding {111} facet. Note that this is equivalent to keep a constant value for 𝑑 𝑔,𝑡 𝑁,𝑇 {111} in equation ( 1), defined at the time 𝑡 𝑁,𝑇 of the grain nucleation in twin relationship. Let us consider two grains, 𝑔 1 and 𝑔 2 , in twin relationship, their common twin plane is denoted 𝐷 {111} and also defined by its cartesian equation. The position of the twin plane 𝐷 {111} remains unchanged. For a given liquid cell 𝜇 located in the neighbourhood of 𝑔 1 and 𝑔 2 , the capture liquid the cell 𝜇 is allowed for only 𝑔 1 or 𝑔 2 , depending on the spatial position of its centre with respect to 𝐷 {111} . If its centre is located beneath 𝐷 {111} , only the capture by grain 𝑔 1 will be authorized. Capture by 𝑔 2 will occur otherwise. The reader will find schematized details of the algorithm in [START_REF] Pineau | Cellular automaton modelling of silicon crystallization with grains in twin relationship[END_REF].

Results and discussion

Comparisons with a benchmark experiment

The 3D CA model described above is applied to simulate an experiment corresponding to observations obtained by X-ray imaging [START_REF] Riberi-Béridot | Strain building and correlation with grain nucleation during silicon growth[END_REF]. This benchmark experiment is fully detailed in [START_REF] Pineau | Cellular automaton modelling of silicon crystallization with grains in twin relationship[END_REF][START_REF] Riberi-Béridot | Strain building and correlation with grain nucleation during silicon growth[END_REF]. Numerical and process parameters as well as material properties are provided in [START_REF] Pineau | Cellular automaton modelling of silicon crystallization with grains in twin relationship[END_REF]. Figure 1 shows the comparison between the simulation and experimental results. (see [START_REF] Pineau | Cellular automaton modelling of silicon crystallization with grains in twin relationship[END_REF] for more details). Twin grain boundaries apart from the one due to Σ3 twin nucleation are the consequences of grain encounters during growth of A-, B-, C-and D-type grains. This result is consistent with the experimental observations. Prior to grain D nucleation, a total of 5 grooves between B-type and C-type grains were formed but did not lead to nucleation of the D-type grain on the {111} plane of a C-type grain.

Quantitative comparisons between simulation and experiment are summarized in figure 2. The area fraction, the number of grains of types A-D and the mean twin band width of grains of type A-C were determined. For a given grain, the latter is defined by the width between two twinning events of a given grain type in a direction perpendicular to the twinning plane. The area of the initial seed is not considered when computing the grain area fraction since we focus on the grains that formed as a result of twinning nucleation events: this is why for band width computations, only edge twin nucleation is considered. One can notice the overall good between simulations and experimental observations. Area fraction of A-type grains obtained in the simulation (41.3%) is very close the area fraction measured in the experiment (39.8%). The number of grains in the computation is also higher (24) than the one in the experiment (22), the mean band width being consequently a little higher, 0.11 mm in the simulation and 0.10 mm in the experiment. For B-type grains, the area fraction is 4.7% lower in the simulation (33.4%) than in the experiment (38.1%). The number of grains is almost identical and the mean band width in the simulation is lower (0.1 mm) compared to the experiment (0.12 mm). Similar observations can be done for C-type grains.

Application to a cube geometry

The model is now applied at a larger scale, i.e. to a cube of 5×5×5 cm 3 . Three simulation cases are proposed: directional growth in an imposed and vertical temperature gradient with (i) no facet formation nor grains in twin relationship nucleation, (ii) facet formation without grains in twin relationship nucleation and (iii) facet formation and grains in twin relationship nucleation. For the sake of simplicity, these simulations are denoted I, II and III, respectively. Simulation I is similar to dendritic grain growth situation since no facet or twin nucleation are modelled in this case. The purpose of such simulations is to demonstrate the influence of the facet formation solely and formation coupled with twins nucleation on the final grain structure when compared to a dendritic one. In all cases, heterogeneous nucleation takes place on the bottom surface of the cube, the crystallographic orientations of the initial equiaxed grains being randomly chosen.

Simulation I.

We firstly disable facet formation. Since twin nucleation can occur only on facets, these specific events are also disabled. The corresponding results are shown in figure 3(a). The final grain structure in displayed on the left part of figure 3(a) with a cold colour related to the misorientation, 𝛼 𝑥 , with respect to the x direction (growth direction). It corresponds to the lowest angle formed between all <100> crystallographic directions of the grains and the x axis. The centred part of figure 3(a) shows the IPF map in the x direction and the right one gives the CSL map. The warm colours in figure 3(a) correspond to highly misoriented grains as shown at the very beginning of the simulation at nucleation stage. The growth of highly misoriented grain is blocked in all <100> crystallographic directions. Cold colours are the only one to remain at the end of the solidification which indicates that well orientated grains, with respect to the temperature gradient, can survive. This is simply due to the grain growth competition which commonly occurs between columnar grains during directional solidification.

Simulation II.

Facet formation is now enabled and nucleation of grains in twin relationship is disabled. The corresponding results are shown in figure 3(b). In comparison with figure 3(a), less grain remains at the end of the solidification. Regarding the misorientation, the colour code indicates that more misoriented grains survived (warm colours). This comes from the fact that some of the grains develop {111} facets and become constrained during their growth. Neighbouring grains, that grow without developing {111} facets, can solidify into the free liquid present in the facet and potentially block growth of side grains that form {111} facets. For simulations I and II, the grain density 𝑛 𝑔 decreases continuously. However, the value at time 𝑡 = 4000 s is lower in simulation II than in simulation I demonstrating that grains development is stopped due to the formation of {111} facets on some of the grains. The value of 𝛼 ̅ 𝑥 at 𝑡 = 4000 s is higher in simulation II when compared to simulation I. This is in agreement with figure 3 where more misoriented grains survive as shown in case II. Simulation III includes nucleation of grains in twin relationship. An increase in grain density at 𝑡 ≈ 2600 s is shown. It indicates that rate of twin grains nucleation overcomes grain elimination induced by grain growth competition.

Conclusion

A 3D CA model dedicated to the modelling of silicon solidification in multi-crystalline structure is developed. IPF and CSL maps are simulated and compared with their counterpart produced on a fully solidified sample from a benchmark experiment. Qualitative and quantitative comparisons show good agreements and similar tendencies in microstructure evolution. Application to a larger scale demonstrates that facet formation and grains in twin relationship drastically modify the final structure compared to dendritic development. Such numerical tool is suitable to follow microstructure evolution in casting process dedicated to silicon material. Future development should consequently permit grain structure optimization by direct simulation at the process scale. Further developments of the model could include a thermomechanical analysis during solidification of the grain structure, thus assessing the additional stress field that may play a role on the formation of the twin grain boundaries and of dislocations.

Figure 1 (

 1 a) is the inverse pole figure (IPF) map along the growth direction and displayed on the longitudinal section (left) simulated and (right) measured.Figure 1(b) shows the coincidence site lattice (CSL) map (left) simulated and (right) measured. The two sets of Σ3 twin relationships with a common (11 ̅ 1) plane on the right-hand side (A/B type grain boundaries) and (111 ̅ ) plane on the left-hand side (A/C type grain boundaries) are made of grains with various band widths. This is directly controlled by the value of the standard deviation to activate the {111} facets, Δ𝑇 𝑔,𝑠 {111}

Figure 1 .

 1 Figure 1. Qualitative comparison between the simulation results and the benchmark experiment: (a) IFP maps along the growth direction on longitudinal sections and (b) CSL maps revealing the different types of twinned GBs.

Figure 2 .

 2 Figure 2. Statistical analyses showing (a) the area fraction, (b) the number of grains and (c) the mean band width between twin boundaries for the simulation results and the benchmark experiment shown in figure 1.

Figure 3 .

 3 Figure 3. Results of simulations corresponding to (a) simulation I with a dendritic growth regime, (b) simulation II with growth restricted to facet formation, (c) simulation III with facet formation and grain nucleation in twin relationship. From left to right: final grain structure observed at time 4000 s showing the misorientation, 𝛼 𝑥 , the centered figure is the IFP map along the growth direction, x, and the CSL map is on the right.

3. 2 . 3 .

 23 Simulation III. Both facet formation and grains in twin relationship nucleation are enabled. The corresponding results are shown in figure3(c). Grains in twin relationship can be clearly observed in final structure. Their presence being also shown in IPF and CSL maps. As it is commonly observed in industrial casting ingots[START_REF] Voigt | Proc 14th Europ. PV Solar Energy Conf[END_REF][START_REF] Voigt | Grain orientation and grain boundaries in cast multicrystalline silicon[END_REF], most of the grains are in Σ3 twin relationship. Similar explanations can be developed (figure4) regarding the evolutions of the grain area density, 𝑛 𝑔 (figure4(a)), and mean misorientation, 𝛼 ̅ 𝑥 (figure 4(b)) with respect to the solidification time, 𝑡.

Figure 4 .

 4 Figure 4. Evolution of (a) the grain area density, 𝑛 𝑔 , and (b) the mean misorientation, 𝛼 ̅ 𝑥 , as a function of the solidification time, 𝑡 for the three simulations cases.

  The coordinates (𝑥 𝜇 , 𝑦 𝜇 , 𝑧 𝜇 ) of the center, 𝐶 𝜇 , of cell  are used in the equation associated to activated facets of grain 𝑔 to compute the quantity 𝑎 𝑔 {111} 𝑥 𝜇 + 𝑏 𝑔 {111} 𝑦 𝜇 + 𝑐 𝑔 {111} 𝑧 𝜇 + 𝑑 𝑔,𝑡 0 {111} . The capture is considered only if a negative value is obtained and inhibited otherwise. This condition corresponds to the facet overpassing the cell center.

2.2. Nucleation of a new grain in twin relationship

Nucleation undercooling of twinned grains on {111} facets is defined at sample edges, Δ𝑇 𝑁,𝑇 𝐸 , and at grain boundary grooves, Δ𝑇 𝑁,𝑇 𝐵 (i.e. in the bulk material) in accordance with experimental observations. Indeed twin formation upon directional solidification is only of Σ3 type and Δ𝑇 𝑁,𝑇 𝐵 has a lower value compared with Δ𝑇 𝑁,𝑇 𝐸

Acknowledgements

The present work was funded by the project CrySaLID (ANR-14-CE05-0046-01) of the French National Research Agency (ANR).