Valérie Geoffroy

Geneviève Payette

Florian Mauffrey

Livie Lestin

Philippe Constant

Richard Villemur
email: richard.villemur@iaf.inrs.ca

Strain-level genetic diversity of Methylophaga nitratireducenticrescens confers plasticity to denitrification capacity in a methylotrophic marine denitrifying biofilm

Keywords: Subjects Biotechnology, Genomics, Marine Biology, Microbiology, Biogeochemistry Denitrification, Methylophaga, Methylotroph, Biofilm, Microbiology, Genomic heterogeneity, Marine environment, Bioprocess

à 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

Denitrification describes the successive reduction of nitrate (NO - 3 ) to nitrite (NO - 2 ), nitric oxide (NO), nitrous oxide (N 2 O), and nitrogen (N 2 ) [START_REF] Van Spanning | The nitrogen cycle: denitrification and its relationship to N 2 fixation[END_REF]. This process is used by bacteria for respiration in environments with low oxygen concentrations and with NO - 3 , NO - 2 , NO, and N 2 O as electron acceptors. The process is driven by metalloenzymes NO - 3 reductase, NO - 2 reductase, NO reductase, and N 2 O reductase [START_REF] Einsle | Structural basis of denitrification[END_REF]. As a facultative trait, denitrification occurs frequently across environments, and is performed by bacteria of diverse origins [START_REF] Zumft | Cell biology and molecular basis of denitrification[END_REF]. Furthermore, numerous microorganisms carry incomplete denitrification pathways for either their growth or detoxification, as NO - 2 and NO are deleterious molecules [START_REF] Kaspar | Nitrite reduction to nitrous-oxide by Propionibacteria-Detoxication Mechanism[END_REF][START_REF] Poole | Nitric oxide and nitrosative stress tolerance in bacteria[END_REF][START_REF] Schreiber | Nitric oxide and nitrous oxide turnover in natural and engineered microbial communities: biological pathways, chemical reactions, and novel technologies[END_REF][START_REF] Simon | Diversity and evolution of bioenergetic systems involved in microbial nitrogen compound transformations[END_REF].

Several studies were carried out in our laboratory on a naturally-occurring multispecies denitrifying biofilm that has developed in a methanol-fed, fluidized denitrification continuous system that treated recirculating water of a marine aquarium at the Montreal Biodome. The biofilm consists of at least 15 bacterial species and numerous protozoans (Labbé et al., 2003;[START_REF] Laurin | Microeukaryote diversity in a marine methanol-fed fluidized denitrification system[END_REF], among which Methylophaga spp. and Hyphomicrobium spp. compose more than 50% of the biofilm [START_REF] Labbé | Microbiological community structure of the biofilm of a methanol-fed, marine denitrification system, and identification of the methanol-utilizing microorganisms[END_REF]. [START_REF] Rissanen | Methylophaga and Hyphomicrobium can be used as target genera in monitoring saline water methanolutilizing denitrification[END_REF] also observed the combination of Methylophaga spp. and Hyphomicrobium spp. in the fluidized-bed type denitrification reactors treating the recirculating seawater of the public fish aquarium SEA LIFE at Helsinki, Finland, suggesting the importance of these two genera in marine denitrification processes.

We isolated Hyphomicrobium nitrativorans NL23 and Methylophaga nitratireducenticrescens JAM1 from the biofilm of the Biodome denitrification system and showed they are the main actors of the denitrifying activities (Labbé et al., 2003;[START_REF] Auclair | Dissimilatory reduction of nitrate in seawater by a Methylophaga strain containing two highly divergent narG sequences[END_REF][START_REF] Auclair | Functional diversity in the denitrifying biofilm of the methanol-fed marine denitrification system at the Montreal Biodome[END_REF]Martineau et al., 2013a;[START_REF] Villeneuve | Methylophaga nitratireducenticrescens sp. nov. and Methylophaga frappieri sp. nov., isolated from the biofilm of the methanol-fed denitrification system treating the seawater at the Montreal Biodome[END_REF]. H. nitrativorans NL23 can carry out complete NO - 3 reduction to N 2 , whereas M. nitratireducenticrescens JAM1 can perform three out of the four reduction steps; it misses the reduction of NO - 2 to NO (Martineau, Mauffrey & Villemur, 2015;Martineau et al., 2013b;[START_REF] Mauffrey | Denitrifying metabolism of the methylotrophic marine bacterium Methylophaga nitratireducenticrescens strain JAM1[END_REF]Mauffrey, Martineau & Villemur, 2015;[START_REF] Villeneuve | Complete genome sequences of Methylophaga sp. strain JAM1 and Methylophaga sp. strain JAM7[END_REF]. In addition to H. nitrativorans NL23 and M. nitratireducenticrescens JAM1, we have demonstrated that a small proportion of other denitrifiers exists in the biofilm of the Biodome denitrification system [START_REF] Auclair | Functional diversity in the denitrifying biofilm of the methanol-fed marine denitrification system at the Montreal Biodome[END_REF], suggesting the potential plasticity of the biofilm in adapting to environmental changes.

Methylophaga species are halophilic marine methylotrophic gammaproteobacteria that do not grow in the absence of NaCl and use one-carbon compounds (e.g., methanol) as sole carbon and energy sources with carbon assimilation proceeding via the 2-keto-3-deoxy-6phosphogluconate (KDPG)-variant ribulose monophosphate (RuMP) pathway [START_REF] Villeneuve | Methylophaga nitratireducenticrescens sp. nov. and Methylophaga frappieri sp. nov., isolated from the biofilm of the methanol-fed denitrification system treating the seawater at the Montreal Biodome[END_REF][START_REF] Boden | Emended description of the genus Methylophaga Janvier, others, 1985[END_REF]. The operating conditions of the Biodome denitrification process led to the enrichment of M. nitratireducenticrescens JAM1 in the biofilm, which is the only reported isolated Methylophaga species that can grow under anoxic conditions with NO - 3 or N 2 O as sole electron acceptors [START_REF] Auclair | Dissimilatory reduction of nitrate in seawater by a Methylophaga strain containing two highly divergent narG sequences[END_REF][START_REF] Mauffrey | Denitrifying metabolism of the methylotrophic marine bacterium Methylophaga nitratireducenticrescens strain JAM1[END_REF]Mauffrey, Martineau & Villemur, 2015). This trait is correlated with the presence in the genome of two gene clusters encoding dissimilatory NO - 3 reductases (narGHJI, referred as nar1 and nar2), two gene clusters encoding cytochrome bc-type complex NO reductases (referred as cnor1 and cnor2) and one gene cluster encoding a dissimilatory N 2 O reductase (nos) [START_REF] Villeneuve | Methylophaga nitratireducenticrescens sp. nov. and Methylophaga frappieri sp. nov., isolated from the biofilm of the methanol-fed denitrification system treating the seawater at the Montreal Biodome[END_REF]Mauffrey, Martineau & Villemur, 2015). M. nitratireducenticrescens JAM1 lacks a gene encoding a dissimilatory copper-(NirK) or cytochrome cd1-type (NirS) NO - 2 reductase.

Our group has also investigated the potential of using the batch operating mode for the denitrification system instead of the continuous one as it prevailed at the Montreal Biodome. With this mode, a better real-time control was obtained to avoid sulfate-reduction that can be easily occurring in seawater biotreatments [START_REF] Labelle | Optimisation d'un bioprocédé de dénitrification d'un sytème aquacole marin en circuit fermé au Biodôme de Montréal[END_REF]. In this report, we aim to assess the impact of the batch operating mode on the bacterial population of the denitrifying biofilm. We observed the disappearance of H. nitrativorans NL23 without the loss of the denitrifying activities during the biofilm acclimation to this mode. We hypothesized that new denitrifiers were enriched during the biofilm acclimation, which displaced H. nitrativorans NL23. From this biofilm, a new denitrifying bacterial strain, named GP59, was isolated that is closely related to M. nitratireducenticrescens JAM1. We compared strain GP59 with M. nitratireducenticrescens JAM1 at the physiological, genomic and transcriptomic levels. Genomic data were also compared with available genomes of Methylophaga species. Our results present new insights of M. nitratireducenticrescens under anoxic environments. They also provide indications that strain-level genomic heterogeneities can influence the dynamism of the biofilm microbial population in adapting to new operating conditions.

MATERIAL AND METHODS

NO -

, NO - 2 and N 2 O measurements

Measurements of NO - 3 , NO - 2 and N 2 O concentrations in all biofilm and planktonic cultures were performed using ion chromatography as described in [START_REF] Mauffrey | Denitrifying metabolism of the methylotrophic marine bacterium Methylophaga nitratireducenticrescens strain JAM1[END_REF].

Biofilm acclimation

Artificial seawater (ASW) medium was composed of (for 1 liter solution): 27.5 g NaCl, 10.68 g MgCl 2 *6H 2 O, 2 g MgSO 4 *7H 2 O, 1 g KCl, 0.5 g CaCl 2 , 456 µL of FeSO 4 *7H 2 O 4 g/L, 5 mL of KH 2 PO 4 51.2 g/L, 5 mL of Na 2 HPO 4 34 g/L. It was supplemented with 1 mL of trace elements (master solution: FeSO 4 *7H 2 O 0.9 g/L, CuSO 4 *5H 2 O 0.03 g/L and MnSO 4 *H 2 O 0.2 g/L) (Labbé, Parent & Villemur, 2003) and NaNO 3 (21.4 mM or 300 mg-N/L). The pH was adjusted (NaOH) at 8.0. Three 120-mL vials containing each 20 unused ''Bioflow 9 mm'' carriers (Rauschert, Steinwiessen, Germany) and 60 mL of ASW were purged of oxygen for 10 min with pure nitrogen gas (N 2 , purity > 99.9%; Praxair, Mississauga, ON, Canada), sealed with sterile septum caps and autoclaved. Prior use, these carriers were washed with HCl 10% (v/v) for 3 h, rinsed with water and autoclaved. A 90 µL volume (72 mg) of filtered-sterilized methanol was then added to the vials (0.15% v/v final concentration).

Carriers (Bioflow 9 mm) with the denitrifying biofilm were taken from the denitrification system at the Montreal Biodome and frozen at -20 • C in seawater with 20% glycerol [START_REF] Laurin | Long-term storage conditions for carriers with denitrifying biomass of the fluidized, methanol-fed denitrification reactor of the Montreal Biodome, and the impact on denitrifying activity and bacterial population[END_REF] until use. The biomass from several frozen carriers was thawed, scrapped, weighted and dispersed in the ASW medium at 0.08 g (wet weight)/mL. The biomass (5 mL/vial, 0.4 g of biofilm) was then distributed with a syringe and an 18 G 1 1 2 needle to three vials. The vials were incubated at room temperature and 100 rpm with an orbital shaker in the dark. In average once a week, the vial was opened in ambient air, the carriers were taken, gently washed with ASW medium to remove the excess medium and the free bacteria, then transferred into vials with fresh anoxic medium, and incubated in the same conditions (Fig. S1). Methanol and NO - 3 were added when needed if NO - 3 was completely depleted during the week. After the 5th transfer, the biofilm was collected. One part was preserved in 15% glycerol at -70 • C for further used such as DNA extraction (see below). The other part was used to isolate bacterial isolates (next subsection).

Isolation of Methylophaga sp. strain GP59

The acclimated biofilm was dispersed in saline solution (3% NaCl, 34.2 mM phosphate buffer pH 7.4), and serial dilutions were made and inoculated onto the Methylophaga medium 1403 (American Type Culture Collection [ATCC], Manassas, VA, USA) supplemented with 1.5% agar and 0.3% v/v methanol. The plates were incubated at 30 • C for no more than 7 days, under oxic conditions (no NO - 3 ) or under anoxic conditions in anoxic jars with a gaspak (BBL TM GasPak TM Plus Anaerobic System). For the latter conditions, the medium was supplemented with 21.4 mM NaNO 3 . Several colonies were picked and cultured in the Methylophaga 1403 liquid medium under oxic or anoxic conditions for 2-7 days. After restreaking three times onto agar medium, isolates were tested for growth under anoxic conditions. NO - 3 and NO - 2 consumption was measured after four days of incubation. Identification of the isolates was determined by extracting total DNA from the isolates, followed by PCR amplification of the 16S ribosomal RNA (rRNA) genes (forward primer 5 -AGAGTTTGATCCTGGCTCAG-3 and reverse primer 5 -AAGGAGGTGATCCAGCCGCA-3 , which correspond to positions 8 to 27 and 1,521 to 1,540 in the Escherichia coli 16S rRNA gene) and the sequencing of the resulting PCR products (McGill University and Genome Quebec Innovation Center, Montréal QC, Canada).

Planktonic pure cultures

M. nitratireducenticrescens GP59 and JAM1 were cultured in the Methylophaga medium 1403 as described by [START_REF] Villeneuve | Methylophaga nitratireducenticrescens sp. nov. and Methylophaga frappieri sp. nov., isolated from the biofilm of the methanol-fed denitrification system treating the seawater at the Montreal Biodome[END_REF] with 0.3% methanol. When required, NO - 3 (NaNO 3 ) was added to the medium. For the anoxic cultures, 70-ml vials with 30 ml medium were flushed with pure nitrogen gas for 20 min, sealed with sterile septum caps and autoclaved. Culture bottles were incubated at 30 • C (unless stated otherwise) in the dark. In the determination of optimal growth conditions, the pH, temperature and NaCl concentrations were adjusted in the Methylophaga medium 1403. For N 2 O measurements, 720-mL bottles with 60-ml medium were used. For oxic cultures, cultures were performed in Erlenmeyer flasks and shaken at 150 rpm and 30 • C [START_REF] Mauffrey | Denitrifying metabolism of the methylotrophic marine bacterium Methylophaga nitratireducenticrescens strain JAM1[END_REF]. Bacterial growth was monitored by spectrophotometry (OD 600 nm ). Bacterial flocs were dispersed with a Potter-Elvehjem homogenizer prior to measurements when needed.

DNA extraction of the biofilm and quantitative PCR

The frozen samples of the biofilm taken from the denitrification system and of the acclimated biofilm taken after the 5th transfer were thawed, and washed 3 times with 0.5 mL the TEN buffer (50 mM Tris-HCl pH 8.0, 10 mM EDTA.2Na, NaCl 150 mM) and 5,000 g centrifugation for 2 min. The biomass was dispersed in 0.5 mL TEN with lysozyme (2.5 mg/mL) and incubated at 37 • C, 30 min. After adding sodium dodecyl sulfate and β-mercaptoethanol (final concentrations of 2% and 1%, respectively), the biomass was subjected to three cycles of freeze-thaw (10 min dry ice/ethanol bath, followed by 5 min at 65 • C), then treated with proteinase K (50 µg/mL final concentration) at 45 • C for 2 h, and finally with RNAse A (20 µg/mL final concentration) at 37 • C for 15 min. The DNA was purified by extractions with phenol/chloroform/isoamyl alcohol (25:24:1 v/v) and chloroform/isoamyl alcohol (49:1), and by precipitation with 0.25 volume of 10 M ammonium acetate and 2 volumes of ethanol 95% (Sambrook & Russell, 2001). DNA concentration was determined by the quantiFluor dsDNA system (Promega, Madison, USA). Quantitative PCR (qPCR) was performed with SYBR green in 20-µL volume with the Fast Start essential DNA Green Master (Roche Diagnostics, Laval, QC, Canada) containing 200 µM of primers (Table 1) and 50 ng biofilm DNA. The amplifications were performed at 95 • C for 10 min, followed by 40 cycles at 95 • C for 20 s, at the annealing temperature (Table 1) for 20 s, and at 72 • C for 15 s. After amplifications, specificity of the PCR products and the presence of primer dimers were verified by performing a melt curve by increasing the temperature from 65 • C to 95 • C by increments of 1 • C per step with a pause of 5 s for each step. Reactions were performed in a Rotor-Gene 6000 real-time rotary analyzer (Qiagen Inc. Toronto, ON, Canada). The amplification efficiencies for all primer pairs varied between 0.9 and 1.1. The copy number of each gene per ng of biofilm was calculated according to standard curves using dilutions of a gel-purified, PCR-amplified fragment of the corresponding genomic region of strain JAM1 (tagH ) or of strain GP59 (nirK ) (Table 1), or using dilution of the NL23 genomic DNA (napA) or the JAM1 genomic DNA (narG1) as the template.

Genome sequencing and annotations

High molecular weight total DNA was extracted from GP59 and JAM1 pure cultures with the DNeasy Blood & Tissue Extraction Kit according to the manufacturer (Qiagen, Hilden, Allemagne). DNA samples were sent to the Genome Quebec Innovation Center for genome sequencing and assembly by using the PacBio technology. Genome assembly used was the HGAP workflow version April 2016 [START_REF] Chin | Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data[END_REF]. Coverage was over 260 times with both genomes. Genome assembly was checked manually with the DNASTAR software (version 14.1) (Madison, WI, USA) for confirming or correcting pseudogenes. Gene annotations for the JAM1 genome were updated accordingly (Genbank accession number CP003390.3). For the GP59 genome, gene annotations were performed at NCBI (GenBank accession number CP021973). Annotations were also performed for both genomes at the Rapid Annotation using Subsystem Technology (RAST) web site (http://rast.nmpdr.org/seedviewer.cgi) and used for further analyses [START_REF] Aziz | The RAST server: rapid annotations using subsystems technology[END_REF][START_REF] Overbeek | The SEED and the Rapid Annotation of microbial genomes using Subsystems Technology (RAST)[END_REF].

Alignment between GP59 and JAM1 genomes was performed at RAST and with the Mauve Multiple Genome Alignment version 20150226 (progressive alignment) [START_REF] Darling | Mauve: multiple alignment of conserved genomic sequence with rearrangements[END_REF]. The resulting alignments were corrected manually to precisely determine the common and unique genes, but also the intergenic sequences. Riboswitches and non coding RNA (ncRNA) sequences were determined by the NCBI annotation system. The average nucleotide identity (ANI) analysis was performed at http://enve-omics.ce.gatech.edu/ani/ with window size 1,000 bp, step size 200 bp, minimum identity 70%, minimum length 700 bp, minimum alignment 50 [START_REF] Goris | DNA-DNA hybridization values and their relationship to whole-genome sequence similarities[END_REF][START_REF] Rodriguez | The enveomics collection: a toolbox for specialized analyses of microbial genomes and metagenomes[END_REF]. Genomes were analysed for tandem repeats at https://tandem.bu.edu/trf/trf.html [START_REF] Benson | Tandem repeats finder: a program to analyze DNA sequences[END_REF].

NirK Phylogenetic analysis

The deduced amino acid sequence of the GP59 nirK was compared with public protein databases by BLASTP at the NCBI web site to retrieve the closest NirK sequences. These sequences were than aligned with Cobalt [START_REF] Papadopoulos | COBALT: constraint-based alignment tool for multiple protein sequences[END_REF]. Evolutionary analysis was conducted using MEGA6.06 [START_REF] Tamura | MEGA6: molecular evolutionary genetics analysis version 6.0[END_REF] with the Maximum Likelihood method based on the Le_Gascuel_2008 model [START_REF] Le | An improved general amino acid replacement matrix[END_REF]. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model, and then selecting the topology with superior log likelihood value. A discrete Gamma distribution was used to model evolutionary rate differences among sites (five categories [+G, parameter = 0.8954]). The analysis involved 21 amino acid sequences. There were a total of 378 positions in the final dataset.

The transcriptomes of M. nitratireducenticrescens GP59 and JAM1

Anoxic cultures of strains GP59 and JAM1 were performed in triplicate in the Methylophaga medium 1403 supplemented with 21.4 mM NaNO 3 . Cells were harvested when NO - 3 was almost reduced in all cultures. Total RNA extraction, RNA processing, cDNA sequencing and transcriptome analyses were described by Mauffrey, Martineau & Villemur (2015), with the exception that we used Transcripts Per Million (TPM) for normalization instead of Reads Per Kilobase Million. One replicate culture of strain GP59 generated inconsistent results and was not considered in the transcriptome analyses. Intergenic regions were also considered in the transcriptome analyses to detect possible missed annotated genes or non-coding expressed sequences. An intergenic region includes the sequence from the stop codon of a gene to the start/stop codon of the adjacent gene. This sequence was then trimmed of 50 nucleotides at both ends to exclude as much as possible the transcript reads that would belong to the adjacent genes. Significance for difference in gene expression of the corresponding genes/sequences (defined as TPM) between strain JAM1 (triplicate cultures) and strain GP59 (duplicate cultures) was performed with the R Bioconductor NOIseq package v2.14.0 (NOIseqBio with replicates and q = 0.9) [START_REF] Tarazona | Differential expression in RNA-seq: a matter of depth[END_REF] and run with the R software v3.2.3 (R Core [START_REF] Core | A language and environment for statistical computing[END_REF]. Genes/sequences of one strain that had >2-fold higher level of TPM than the other strain showed significant differences. RNA-seq reads were deposited in the Sequence Read Archive (SRA) under the numbers SRP066381 (strain JAM1) and SRP132510 (strain GP59) at NCBI.

RESULTS AND DISCUSSION

Biofilm acclimation to batch cultures

The biofilm of the denitrification system (operated under continuous mode) was dispersed in vials containing 20 Bioflow carriers, the same used in the denitrification system, and cultured under anoxic, batch-mode conditions with NO - 3 . Only the carriers with the The biofilm was batch-cultured under anoxic conditions. After the 5th transfer, the NO - 3 and NO - 2 concentrations were measured over different time intervals. NO x : combined measurements of NO - 3 and NO - 2 . Each point is the average with standard deviation of triplicate cultures.

Full-size DOI: 10.7717/peerj.4679/fig- 1 growing biofilm were then transferred each week in fresh medium and cultured under the same conditions (see 'Material & Methods' Fig. S1). At the end of the 5th transfer culture when the NO - 3 and NO - 2 consumption rates stabilized, the buildup of biofilm was enough (Fig. S2) to extract total DNA and to perform culture assays on agar plates.

NO - 3 and NO - 2 concentrations were measured at regular interval times during the 5th transfer culture. Complete reduction of NO - 3 and NO - 2 (NO x ) was achieved within 15 h, at 1.8 mM NO x h -1 (Fig. 1). This is 5-times faster than what was measured by Laurin et al. (24) with 60 Bioflow carriers taken directly from the denitrification reactor that operated under continuous-mode conditions (0.37 mM NO x h -1 or 5.2 NO x -N L -1 h -1 ). This result suggests that changes occurred in the microbial population of the batch biofilm cultures. qPCR assays were performed to determine the concentrations of M. nitratireducenticrescens JAM1 (narG1-targeted primers) and H. nitrativorans NL23 (napA-targeted primers) in the original and acclimated biofilm. In the original biofilm, strain JAM1 was at 2.5 ± 0.6 ×10 4 narG1-copies/ng biofilm DNA and H. nitrativorans NL23 at 6.2 ± 0.4 ×10 5 napA-copies/ng biofilm DNA. In the acclimated biofilm, the concentration of strain JAM1 increased by one order of magnitude (3.8 ± 0.9 ×10 5 narG1-copies/ng biofilm DNA), but decreased substantially for strain NL23 (2.0 ± 0.5 × 10 2 napA-copies/ng biofilm DNA). With the decrease of strain NL23 and the increase of strain JAM1 in the acclimated biofilm with a gain of denitrification performance, we hypothesized that a subpopulation of M. nitratireducenticrescens that can perform the complete denitrification pathway was enriched in the acclimated biofilm.

Strain GP59 displaced H. nitrativorans NL23 and M. nitratireducenticrescens JAM1

The acclimated biofilm from the 5th transfer cultures was dispersed and cultured onto a Methylophaga specific medium, and several isolates were screened for NO - 3 and NO - 2 consumption (Fig. S1). One isolate, which we named strain GP59, was capable to consume NO - 3 and NO - 2 , and cultures showed gas production, suggesting full denitrifying activities by this strain. The 16S rRNA gene sequence of strain GP59 has 100% identity to that of M. nitratireducenticrescens JAM1.

The concentration of strain GP59 was measured in the acclimated biofilm by qPCR to determine whether this strain is a major denitrifier. To distinguish strain GP59 from strain JAM1, primers targeted strain-specific sequences (tagH for strain JAM1 and nirK for strain GP59) were developed based on their respective genome (Table 1; see below). The GP59 cell concentrations went from undetected (below limit of PCR detection) in the original biofilm to high level in the acclimated biofilm at 2.2 ± 0.35 ×10 5 nirK -copies/ng biofilm DNA. This is 2-3 orders of magnitude higher than that of strain JAM1 (4.0 ± 1.7 ×10 2 tagH -copies/ng biofilm DNA). These results confirm our hypothesis that a new M. nitratireducenticrescens strain with full denitrification capacity occurred in the acclimated biofilm, and displaced H. nitrativorans NL23, but also M. nitratireducenticrescens JAM1. Comparative analyses were then performed between strain GP59 and M. nitratireducenticrescens JAM1 at the physiological, genomic and transcriptomic levels.

Physiological characterization of strain GP59

Strain GP59 cultured under anoxic conditions (planktonic pure cultures) showed a 24-h lag before growth occurred (Fig. 2A). The growth yields reached about 1.2 OD 600nm with 42.8 mM NO - 3 exposure and did not increase significantly in cultures exposed with higher NO - 3 concentrations. Strain JAM1 cultured under the same conditions showed no lag phase, and growth yields were 3-to 6-fold lower than those of the GP59 cultures (Fig. 2A). The maximum specific growth rate (µmax) and the half-saturation constants of NO - 3 for growth (Ks) were calculated (Table 2). The µmax for strain GP59 is higher (3.3 times) than that of strain JAM1. To assess the affinity of strain GP59 toward NO - 3 for growth, the µmax/Ks ratio was calculated [START_REF] Healey | Slope of the monod equation as an indicator of advantage in nutrient competition[END_REF] (Table 2). This ratio (1.2 µM -1 NO - 3 h -1 ) is not different than the one calculated for strain JAM1 at 1. Martineau & Villemur, 2015). These results concur with the genome sequences (see below) with near 100% identity between the two strains in gene clusters encoding the two Nar systems and the NO - 3 transporters (NarK). As observed with growth, a 24-h lag period was observed in the GP59 cultures before NO - 3 started to be consumed, whereas complete NO - 3 consumption occurred in the JAM1 cultures within 24 h (Fig. 2B). The NO - 3 reduction rates increased linearly with the increase of NO - 3 concentrations in the GP59 cultures (Fig. 3A). The specific NO - 3 reduction rates (rates normalized by the biomass) averaged around 1.5 to 2 mM NO - 3 h -1 OD -1 and showed no significant changes at any NO - 3 concentrations tested (Fig. 3B). In the JAM1 cultures, the NO - 3 reduction rates reached a plateau at 24 mM NO - 3 (Fig. 3A), and, as observed with the GP59 cultures, showed no significant changes in the specific NO - reduction rates (averaged around 4 to 6 mM NO - 3 h -1 OD -1 ) (Fig. 3B). Interestingly, these specific rates are 2.4 to 5.4-fold higher than those of strain GP59, suggesting that the JAM1 cells have a higher dynamism of NO - 3 processing (e.g., NO - 3 intake and reduction) than the GP59 cells. Close examination of the respective genomes and transcriptomes (see below) did not reveal specific gene(s) that would explain these differences.

3 µM -1 NO - 3 h -1 (Mauffrey,
N 2 O was detected in the GP59 anoxic cultures and reached maximum accumulation (ca. 0.5 µmol N 2 O vial -1 or 0.04% N-input) when NO - 2 peaked in the medium, and then decreased in concentration afterwards (Fig. 2B). This result concurs with the presence of a µmax/Ks (µM

-1 h -1 ) 1.2 (0.2) 1.3 (0.3) Notes.
µmax, maximum growth rates. Ks, half-saturation constants of NO - 3 for growth. Values between parentheses are standard deviation of triplicates. a From Mauffrey, Martineau & Villemur (2015). 3 reduction rates. These rates were calculated with the NO - 3 reduction rates by the generated culture biomass (OD 600nm ) at the end of the exponential phase.

Full-size DOI: 10.7717/peerj.4679/fig- 3 complete denitrification pathway in strain GP59 as revealed by its genome sequence (see below). Strain GP59 was cultured under oxic conditions with different concentrations of NO - 3 (9.6 to 142.8 mM). Contrary to the JAM1 cultures that can reduce NO - 3 to NO - 2 under oxic conditions (Mauffrey, Martineau & Villemur, 2015), no NO - 3 consumption was observed in the GP59 cultures at any tested concentrations. These latter results, coupled with the 24-h lag period before nitrate reduction occurs in the GP59 anoxic cultures, suggest that strains GP59 and JAM1 regulate differently the expression of their denitrification genes. Both strains cannot grow on methylamine and fructose. They showed similar growth profile related to NaCl concentrations (optimal growth at 1 to 5% NaCl; weaker growth at 8%; no growth at 0 and 10%), pH (optimal growth at pH 8; weaker growth at pH 7, 9 and 10; no growth at pH 6 and 11) and temperature (optimal growth at 30 • C; variability in culture replicates at 22 and 34 • C; weak growth at 37 • C). Sensibility to antibiotics by strain GP59 was also similar to those reported by [START_REF] Villeneuve | Methylophaga nitratireducenticrescens sp. nov. and Methylophaga frappieri sp. nov., isolated from the biofilm of the methanol-fed denitrification system treating the seawater at the Montreal Biodome[END_REF] for strain JAM1 (chloramphenicol, trimethoprim, ampicillin, tetracycline, streptomycin, gentamycin and kanamycin).

The genome of strain GP59 contains nirK

The genome of strain GP59 was sequenced by the PacBio technology. Although strain JAM1 was already sequenced by the pyrosequencing technology, it was resequenced by the PacBio technology to correct some discrepancies in the genome assembly (Data S1). Analysis of the GP59 genome revealed the presence of two prophages (named GPMu1 and GPMu2) integrated side by side in the same orientation and separated by 11,102 nt (Data S2). The intercalated sequence has a GC content of 40.1%, which is lower than in the overall genome (44.8%). Gene arrangement of the prophages resembles to that of the E. coli phage Mu and the Haemophilus influenzae FluMu prophage [START_REF] Morgan | Bacteriophage Mu genome sequence: analysis and comparison with Mu-like prophages in Haemophilus, Neisseria and Deinococcus[END_REF]. Gene clusters encoding Mu-type prophage proteins are also present in the genome of two Methylophaga species such as M. frappieri JAM7 that was also isolated from the original biofilm of the denitrification system [START_REF] Auclair | Dissimilatory reduction of nitrate in seawater by a Methylophaga strain containing two highly divergent narG sequences[END_REF] (Data S2). No prophage was found in the JAM1 genome.

The 11,102 nt intercalated sequence contains a gene encoding a NirK-type NO - 2 reductase of 363 amino acid residues (Data S2). Upstream of this nirK, is a gene encoding the NO reductase activation protein NorD. The precise function of NorD is yet to be known, but the norD mutant in Paracoccus denitrificans showed no activity of the cd 1 -type NO - 2 reductase (NirS) and of the NO reductase [START_REF] Deboer | Mutational analysis of the nor gene cluster which encodes nitric-oxide reductase from Paracoccus denitrificans[END_REF]. Available genomic sequences from other Methylophaga species and from metagenomic studies were screened for the presence of genes encoding NirK highly similar to the GP59 NirK. The closest affiliations of the GP59 NirK are with NirK encoded in three reconstituted Methylophaga sp. genomes from a metagenomic study of marine subsurface aquifer samples (NCBI Bioproject PRJNA391950 by Tully, Wheat, Glazer and Huber, University of Southern California) (Fig. 4). These three NirK have 84-85% similarity with the GP59 NirK sequence. In addition, norD is also present adjacent to nirK in these three Methylophaga sp. genomes.

The previous results suggest that there are Methylophaga strains in natural environments that carry denitrification modules. We further performed data mining to retrieve denitrification genes in available genomes of Methylophaga species or in reconstituted Methylophaga genomes from metagenomic studies. Data provided by the Bioproject PRJNA391950 revealed that the metagenome of the Methylophaga sp. isolate NORP53 (accession number NVVW00000000) contains two gene clusters with one (contig NVVW01000003) having the same gene arrangement of the Nar1 gene cluster in the An unrooted phylogenetic tree demonstrating the evolutionary relationship of nirK deduced amino acid sequences is illustrated. The evolutionary history was inferred by using the Maximum Likelihood and the tree with the highest log likelihood (-7,598.4857) is shown. The percentage (values above 70%) of trees in which the associated taxa clustered together is shown next to the branches (bootstrap analysis with 1,000 replicates). The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The GenBank accession number is provided after each bacterial name. *: The Cycloclasticus sp. NORP92, Gammaproteobacterium NORP60, and the Methylophaga sp. NORP93, NORP33 and NORP65 are from the NCBI bioproject PRJNA391950 by Tully, Wheat, Glazer and Huber, University of Southern California.

Full-size DOI: 10.7717/peerj.4679/fig- 4 JAM1 and GP59 genomes (narXL, NarK1K2GHJI ), and another one (NVVW01000010) encoding a N 2 O reductase (nosRZDFYL). At the amino acid sequence level, these gene clusters of isolate NORP53 are 67 to 90% similar with the corresponding gene products in the JAM1 and GP59 genomes. Finally, nirK sequences were found in the genome of M. frappieri JAM7 and of 34 Methylophaga sp. retrieved from metagenomic studies (Data S3). Although M. frappieri JAM7 does not grow under anoxic conditions with NO - 3 , which concurs with the absence gene cluster encoding a dissimilatory NO - 3 reductase, the occurrence of nirK may be associated with the detoxification mechanism. All these results suggest that the different denitrification modules are carried by some Methylophaga species present in natural environments. region containing the denitrification gene clusters nar1, nar2, nor1, nor2 and nos (Mauffrey, Martineau & Villemur, 2015) is 99.94% identical (43 substitutions, no gap) between the two strains. Based on the genome annotations, functions associated with CDS that are unique to strain JAM1 have in most cases an equivalent in strain GP59. In addition to CDS associated with the two prophages and the intercalated sequence, the GP59 genome has CDS associated with an integrated plasmid (position 724937 to 765476; Data S4). The GC skew between the two genomes is relatively similar (Fig. S3). Two-way ANI analysis showed 99.30% identity, and >80% conserved sequences between the two genomes, indicating that both strains belong to the same species [START_REF] Goris | DNA-DNA hybridization values and their relationship to whole-genome sequence similarities[END_REF][START_REF] Rodriguez | The enveomics collection: a toolbox for specialized analyses of microbial genomes and metagenomes[END_REF].

Strain GP59 belongs to the species Methylophaga nitratireducenticrescens

Both genomes have small tandem repeats with >10 repeats positioned at equivalent locations. These tandem repeats vary in number of repeats between the two strains (Table 3), and are located in CDS with no putative known function. A long tandem repeat (four 3095-nt repeats) was found in the GP59 genome (Table 3). Genes encoding an aspartate aminotransferase and a TonB-dependent siderophore receptor are among the three genes repeated 4 times. The equivalent region in the JAM1 genome shows no repetition, with 93% identity with the 3095-nt repeat sequence with most of the divergence in the TonB-dependent siderophore receptor gene.

Strain GP59 contains two plasmids here named pGP32 and pGP34 (Table 3). This result concurs with the observation of extra-chromosomic DNA when plasmid DNA extraction was performed on GP59 biomass (Data S5). These plasmids encode type IV secretion system as observed in many plasmids hosted by Gram negative bacteria (Data S5) [START_REF] Christie | Type IV secretion: intercellular transfer of macromolecules by systems ancestrally related to conjugation machines[END_REF]. They also carry genes encoding for recombinases, nucleases, DNA methyltransferases, toxin/antitoxin, and carbon storage regulator. These two plasmids show no homology with the integrated plasmid in the GP59 genome. M. frappieri JAM7 also carries a 48.5 kb plasmid, which resembles in its structure to the IncPalpha plasmids (Data S5). It shares no homology with the pGP32 and pGP34 plasmids.

M. nitratireducenticrescens GP59 did not evolve from M. nitratireducenticrescens JAM1

Did strain GP59 evolve during the acclimation process of the biofilm from strain JAM1 that underwent series of chromosomic rearrangements and phage infections, or was strain GP59 present in the original biofilm and enriched during the acclimation? To try to answer this question, we analyzed the Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR) region that has been identified in the GP59 and JAM1 genomes. The CRISPR system is an adaptive immunity system that is present in most archaea and many bacteria and that acts against invading genetic elements, such as bacteriophages and plasmids. The chromosomic arrangement of CRISPR system includes CRISPR associated genes (cas), a leader sequence upstream of an array of short repeats interspersed with unique spacers that are almost identical to fragments of bacteriophage and plasmid genes. The transcription of the spacer array provides complementary DNA targeting, for instance, phage genome resulting in target DNA degradation [START_REF] Barrangou | CRISPR-Cas systems: prokaryotes upgrade to adaptive immunity[END_REF][START_REF] Garneau | The CRISPR/Cas bacterial immune system cleaves bacteriophage and plasmid DNA[END_REF][START_REF] Makarova | Evolution and classification of the CRISPR-Cas systems[END_REF]. Addition of new spacers in response, for instance, of bacteriophage infection occurs proximal to the leader sequence. Accordingly, CRISPR spacers provide a historical perspective of phage exposure, with spacers at the vicinity of the leader were relatively recently added, and those distal spacers likely originated from previous events [START_REF] Horvath | Diversity, activity, and evolution of CRISPR loci in Streptococcus thermophilus[END_REF].

The CRISPR region in GP59 and JAM1 genomes comprises 5 associated CRISPR genes (Fig. 5A). Except for the spacer sequences, the nucleotide sequence of the CRISPR genes and the repeat sequences are identical in both genomes (Figs. 5A,5C). There are 105 and 115 spacers in the JAM1 and GP59 CRISPR regions, respectively (Fig. 5B). In the JAM1 CRISPR region, 14 spacer sequences were found twice. The GP59 and JAM1 CRISPR regions share 29 spacer sequences (Fig. 5B); all of them are distal to the leader sequence (the most ancient acquired sequences, [START_REF] Horvath | Diversity, activity, and evolution of CRISPR loci in Streptococcus thermophilus[END_REF]. No spacer sequence was found twice in the GP59 CRISPR region. As both CRISPR regions do not share sequence between their proximal spacers, strains GP59 and JAM1 underwent a different history of phage infections from a common ancestry. This strongly suggests that strain GP59 did not originate from strain JAM1 that underwent series of chromosomic rearrangements and phage infections in a short period of time and in a closed environment, but rather was present in the original biofilm and was enriched during the acclimation phase of the biofilm. c d e f g h i j k l m n o p q r s t u v w x y z aa 1 2 3 4 5 6 7 8 9 101112 c d e f g h i j k l m n o p q r s t u v w x y z aa 1 2 3 4 5 6 7 8 9 10 1112 The GP59 and JAM1 CRISPR regions were also compared to other available Methylophaga genomes. M. frappieri JAM7 contains 68 repeat units. The repeat sequence has one-nt substitution with the GP59/JAM1 one (Fig. 5C). None of its spacer sequences were found in the GP59 and JAM1 CRISPR regions. In the M. aminisulfidivorans MP, M. sulfidovorans DSM 11578, M. lonarensis MPL genomes, 59, 6 and 63 spacer sequences were found, respectively; none were repeated twice or were found in other Methylophaga species.

The CRISPR genes and gene arrangement in M. frappieri JAM7 are typical of the subtype I-F CRISPR-Cas system, whereas those of M. aminisulfidivorans MP, M. sulfidovorans DSM 11578, M. lonarensis MPL are related to the subtype I-E one (Fig. 5A) [START_REF] Makarova | Evolution and classification of the CRISPR-Cas systems[END_REF]. In strains GP59 and JAM1, the CRISPR genes and gene arrangement cannot be associated to a defined subtype. However, the gene arrangement and the deduced amino acid sequences of cas3 and of the two downstream CDS are highly similar to CRISPR genes found in several bacterial genomes such as Legionella brunensis (GenBank accession number LNXV01000004) and Moraxella atlantae (GenBank accession number LZMZ01000013).

Transcriptome analyses

The GP59 and JAM1 genomes share more than 85% of their genes. However, gene variations, the presence of strain-specific genes, and the occurrence of plasmids in strain GP59 could affect the expression pattern of these common genes. To assess whether these genes were expressed differently between strains GP59 and JAM1, the transcriptome of both strains was derived from anoxic cultures at the end of the exponential growth when NO - 3 was nearly completely reduced in both cultures. Among the 2857 common genes (CDS and ncRNA), more than 85% of those in strain GP59 showed no significant difference in the transcription levels (<2-fold difference in TPM) (Data S4) with their counterparts in strain JAM1. The vast majority of genes involved in the nitrogen pathways were expressed in equivalent levels in both strains (Table 4). The fact that the cultures were sampled for RNA extraction when both strains were at an equivalent physiological phase (anoxic conditions with NO - 3 almost consumed) could explain this high proportion of similar gene expression between the two strains. Among genes with significant differences are those encoding the synthesis of the osmoprotectant ectoine, which were 3-to 6-fold higher in the JAM1 cultures than in the GP59 cultures (Fig. S4). This suggests that strain JAM1 could be more resilient to changes in salinity.

The most expressed genes were sorted out from both transcriptomes to assess which specific metabolic pathways are the most solicited under denitrifying conditions. High levels of gene expression were found in both strains for genes involved in growth, energy and respiration (e.g., ribosomal proteins, ATP synthase, pseudoazurin), and the C-1 metabolic carbon such as methanol dehydrogenase, ribulose monophosphate enzymes (D-arabino-3-hexulose 6-phosphate formaldehyde-lyase, 6-phospho-3-hexuloisomerase), formaldehyde activating enzymes (Data S6, Fig. S4). One of the highest expressed genes in both genomes is located in the chromosomic region comprising three adjacent genes encoding similar proteins (>88% similarity). Each of these genes has a cyclic di-GMP riboswitch in their respective upstream sequence. The first of the three genes showed high level of TPM in the riboswitch and in the coding sequence (Table 5). The other two genes and their respective riboswitch showed much lower expression levels. The deduced amino acid sequence of these three similar proteins contains a carbohydrate-binding module 6 that is found in carbohydrate-active enzymes such as glycoside hydrolases [START_REF] Henshaw | The family 6 carbohydrate binding module CmCBM6-2 contains two ligand-binding sites with distinct specificities[END_REF], and three transmembrane domains. Upstream and downstream of these three genes are genes encoding proteins involved in polysaccharide export and in the capsular polysaccharide biosynthesis. These proteins may participate in this synthesis. Finally, the pGP32 plasmid accounted for 23.1% of total TPM in the GP59 transcriptomes, whereas it was 0.7% for the pGP34 plasmid (Data S6). In the pGP32 plasmid, nine genes accounted for 20% of total TPM, among which are genes involved in restriction-modification systems (Data S6).

CONCLUSIONS

Culturing the denitrifying biofilm under batch conditions has favored the enrichment of a new denitrifying subpopulation, representing by M. nitratireducenticrescens GP59 that displaced H. nitrativorans NL23 and M. nitratireducenticrescens JAM1. In addition of complete denitrification by strain GP59, two important differences were observed between the two strains. Strain GP59 cannot reduce NO - 3 under oxic conditions as does strain JAM1, and a 24-h lag period was observed in GP59 anoxic cultures before growth and NO - 3 reduction occurs. These results suggest differences in regulation of the denitrification genes between the two strains. The genome of strain GP59 showed insertions of two Mu-type prophages that probably brought a nirK gene, the missing denitrification gene in strain JAM1. Both genomes share more than 85% of their genes, and these genes were expressed at similar level under anoxic conditions at the end of the exponential growth. Finally, analysis of the CRISPR region suggests that strain GP59 did not originate from strain JAM1, but was present in the original biofilm and was enriched during the acclimation process.

These comparative analyses between strains GP59 and JAM1 unveiled that M. nitratireducenticrescens species enriched in the biofilm of the bioreactor encompasses a mosaic population structure. Although plasticity of the genomic landscape of bacterial species resulting to gene rearrangements, point mutations and lateral gene acquisition was observed in nature [START_REF] Allen | Genome dynamics in a natural archaeal population[END_REF]Bendall et al., 2016), the observation of M. nitratireducenticrescens GP59 reported here shows functional consequences of population variation on process rate and system performances. The reason why strain GP59 did occur in the acclimated biofilm but not in the Biodome denitrification system remains obscure and is under investigation. This could be related to the batch-mode conditions where the transient accumulation of NO - 2 (up to 10 mM; Fig. 1) could have adverse the growth of H. nitrativorans NL23, favoring the enrichment of M. nitratireducenticrescens GP59. In the denitrification reactor at the Biodome, such accumulation of NO - 2 was not encountered probably due the continuous-operating mode or the lower denitrification rates that prevailed.

Our study broadens the ecology of M. nitratireducenticrescens with the occurrence of a microbial seedbank suited for anoxic conditions. So far among the genus Methylophaga, M. nitratireducenticrescens is the only reported species that can grow under anoxic conditions. Based on metagenomic studies, however, where denitrification genes were found in Methylophaga-reconstituted genomes, this anoxic metabolism within this genus is more widespread than previously thought. Methanol is generated by anthropogenic activities or by biogenic metabolisms such as plant growth and decomposition, and ocean ecosystems are thought to participate actively in its production. Furthermore, oceans are known to be a sink of atmospheric methanol [START_REF] Galbally | The production of methanol by flowering plants and the global cycle of methanol[END_REF][START_REF] Millet | New constraints on terrestrial and oceanic sources of atmospheric methanol[END_REF]. As they are ubiquitous in marine environments, Methylophaga spp. are therefore important actors in the methanol turnover on Earth, either under oxic environments or, based on our studies with M. nitratireducenticrescens, under anoxic marine environments.

• Livie Lestin performed the experiments, analyzed the data, approved the final draft.

• Philippe Constant conceived and designed the experiments, analyzed the data, contributed reagents/materials/analysis tools, authored or reviewed drafts of the paper.

• Richard Villemur conceived and designed the experiments, performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, prepared figures and/or tables, authored or reviewed drafts of the paper, approved the final draft.

DNA Deposition

The following information was supplied regarding the deposition of DNA sequences:

The genome of strain JAM1 is accessible via GenBank accession number CP003390.3. The genome of strain GP59 is accessible via GenBank accession number CP021973. The RNA-seq reads were deposited in the Sequence Read Archive (SRA) under the numbers SRP066381 (strain JAM1) and SRP132510 (strain GP59) in GenBank.

Figure 1

 1 Figure1NO - 3 and NO - 2 reduction by the acclimated biofilm. The biofilm was batch-cultured under anoxic conditions. After the 5th transfer, the NO - 3 and NO - 2 concentrations were measured over different time intervals. NO x : combined measurements of NO - 3 and NO - 2 . Each point is the average with standard deviation of triplicate cultures.Full-size DOI: 10.7717/peerj.4679/fig-1

Figure 2

 2 Figure 2 Growth, and NO - 3 , NO - 2 and N 2 O dynamics by Methylophaga nitratireducenticrescens GP59 and JAM1. (A) Growth. Strain GP59 and strain JAM1 were cultured under anoxic conditions with 9.6 or 42.8 mM NO - 3 . (B) NO - 3 , NO - 2 and N 2 O dynamics. Strain GP59 and strain JAM1 were cultured under anoxic conditions with 21.4 mM NO - 3 , and NO - 3 , NO - 2 and N 2 O were measured over different time intervals. Each point is the average with standard deviation of triplicate cultures for (A), and of duplicate cultures for (B). Full-size DOI: 10.7717/peerj.4679/fig-2

Figure 3

 3 Figure 3 Specific NO -3 reduction rates. (A) NO - 3 reduction rates of strains GP59 and JAM1. Each point is the average with standard deviation of triplicate anoxic cultures. (B) Specific NO - 3 reduction rates. These rates were calculated with the NO - 3 reduction rates by the generated culture biomass (OD 600nm ) at the end of the exponential phase.Full-size DOI: 10.7717/peerj.4679/fig-3

Figure 4

 4 Figure4Phylogenetic analysis of NirK of Methylophaga nitratireducenticrescens GP59. An unrooted phylogenetic tree demonstrating the evolutionary relationship of nirK deduced amino acid sequences is illustrated. The evolutionary history was inferred by using the Maximum Likelihood and the tree with the highest log likelihood (-7,598.4857) is shown. The percentage (values above 70%) of trees in which the associated taxa clustered together is shown next to the branches (bootstrap analysis with 1,000 replicates). The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The GenBank accession number is provided after each bacterial name. *: The Cycloclasticus sp. NORP92, Gammaproteobacterium NORP60, and the Methylophaga sp. NORP93, NORP33 and NORP65 are from the NCBI bioproject PRJNA391950 by Tully, Wheat, Glazer and Huber, University of Southern California.Full-size DOI: 10.7717/peerj.4679/fig-4

Figure 5

 5 Figure 5 CRISPR chromosomic arrangement in Methylophaga genomes. (A) Gene arrangement. Sequences were retrieved from GenBank. Accession numbers: M. nitratireducenticrescens JAM1: CP003390 (Q7A_2613 to Q7a_2609) and GP59: CP021973 (CDW43_11545 to CDW43_11525); M. frappieri JAM7 CP003380 (Q7C_1045 to Q7C_1040); M. aminisulfidivorans MP: NZ_AFIG01000001 (MAMP_RS08465 to MAMP_RS08425), M. sulfidovorans DSM 11578: FOSH01000003 (SAMN04488079_1032 to SAMN04488079_10310); M. lonarensis MPL, APHR01000054 (MPL1_10022 to MPL1_10062). Nomenclature is based on Barrangou & Marraffini (2014). HP, Hypothetical protein. Numbers at the right are the number of repeat units with the repeat/spacer nucleotide length in parentheses. (B) Spacer/repeat unit arrangement. Blue: the 29 common units. Green: the 14 units repeated twice in strain JAM1. (C) CRISPR repeat sequences in different Methylophaga species. Full-size DOI: 10.7717/peerj.4679/fig-5

Table 1 PCR primers used for qPCR assays. Gene target Sequence (5 -3 ) Annealing temperature ( • C) Fragment length (nt) Comment

 1

	narG1 (strain JAM1/GP59; locus Q7A_446). Nitrate reductase subunit alpha (Nar1 system)
	narG1	AGCCCACATCGTATCAAGCA	61	149	qPCR
	narG1	CCACGCACCGCAGTATATTG
	tagH (strain JAM1, locus tag: Q7A_1110). Teichoic acid export ATP-binding protein a
	TagH Forward	CCGTCATTTCGCTTCAAGAT	55	711	Standard
	TagH Reverse	TCATGGCTTTTTCAGCCTTT
	qTagH Forward	GTTGCAAGGCTATAGTCGGAGT	55	119	qPCR
	qTagH Reverse	TGGTACGCATTCCAGATGAATA
	nirK (strain GP59, locus tag: CDW43_15165). Copper containing nitrite reductase
	NirK Forward	CGTTCAATACATGGGGTAAAGG	55	1,134	Standard
	NirK Reverse	TGGGGCACAGTGATAAACAA
	qNirK Forward	AAGTCGGTAAAGTAGCCGTTGA	55	138	qPCR
	qNirK Reverse	TCTCCATCGTCATTTGAACAAC
	napA (strain NL23, locus tag W911_13875) periplasmic nitrate reductase
	qnapA Forward	AGGACGGGCGGATCAATTTT	61	131	qPCR
	qnapA Reverse	CGGATATGCATCGGACACGA
	Notes.
	a The function was deduced by RAST annotations.
	GenBank accession number of the genome of strain JAM1: CP003390.3, strain GP59: CP021973.1, and strain NL23: CP006912.1.

Table 2 Kinetics of growth under anoxic conditions.

 2

	Strain GP59	Strain JAM1 a

Geoffroy et al. (2018), PeerJ, DOI 10.7717/peerj.4679 12/25

	Nitrosomonas cryotolerans SFP39808.1
	Nitrosomonas mobilis SCZ86535.1
	Nitrosococcus oceani ATCC 19707 ABA56624.1
	Cycloclasticus sp. NORP92 PHS70052.1
	Bathymodiolus azoricus thioautotrophic gill
	symbiont WP_090717277.1
	Methylophaga nitratireducenticrescens GP59, AUZ86105
	Gammaproteobacterium NORP60 PCJ30753.1
	Methylophaga sp. NORP93 PHS69590.1
	Methylophaga sp. NORP33 PHR81400.1
	Methylophaga sp. NORP65 PHS29800.1
	Sunxiuqinia dokdonensis KOH43527.1
	Nitrosovibrio tenuis SEL66906.1
	Nitrosospira multiformis SDQ31436.1
	Nitrosomonas marina SES75436.1
	Shewanella denitrificans OS217 ABE56757.1
	Marinobacter psychrophilus AKO51897.1
	Oleispira antarctica RB-8 CCK74932.1
	Pseudomonas nitroreducens OWP48825.1
	Pseudomonas mosselii SJ10 AIN60781.1
	Rubellimicrobium mesophilum DSM 19309
	EYD77238.1
	Roseovarius pacificus SHM21003.1

  Table 3 summaries the overall features of the GP59 and JAM1 genomes. Both genomes share 2,790 coding sequences (CDS) and 11 riboswitches (Data S4). The 67 kbp chromosomic

Table 3 Genomic features of M. nitratireducenticrescens strains JAM1 and GP59.

 3

	JAM1	GP59

a Based on RAST annotations. b The possible CDS in the ambiguous region of GPMu1 were not taking into account. Common CDS: see Data S4. c Number of copies. Genome location under parenthesis. Only repeats with more than 10 copies are presented. d Number of copies. Genome location under parenthesis. % identity between the 3095-nt long repeats (see text).

Table 4 Transcript levels of genes involved in the nitrogen pathway in Methylophaga nitratireducenticrescens JAM1 and GP59.

 4

	Locus Tag	Gene/function	a Expression ratio		b mxaI = 100
			JAM1/GP59	GP59	JAM1
	Nitric oxide reductase 1
	CDW43_RS01925-40	norCBDQ	0.93	11.7	12.4
	CDW43_RS01945	norR	1.42	4.9	6.8
	CDW43_RS01950	norE	0.84	7.8	6.5
	Nitrate reductase 1 (with transporters and regulators)
	CDW43_RS01975-80	narXL	2.07	3.8	7.9
	CDW43_RS01985-2010	narK1K2narGHJI	0.66	17.5	9.8
	Nitrous oxide reductase
	CDW43_RS02060-85	nosRZDFYL	1.51	5.5	9.0
	Nitrate reductase 2 (with transporter)
	CDW43_RS02155	narK	0.75	1.9	1.4
	CDW43_RS02165-80	narGHJI	0.83	1.9	1.5
	Nitric oxide reductase 2
	CDW43_RS02185-2200	norCBQD	0.83	2.1	1.7
	Copper containing nitrite reductase
	CDW43_RS15160	nirK	na	1.0	na
	Nitrate assimilatory pathway
	CDW43_RS11575	Nitrate/nitrite transporter	1.17	0.2	0.2
	CDW43_RS11585	Assimilatory nitrate reductase large sub.	1.17	0.3	0.4
	CDW43_RS11590	Nitrite reductase [NAD(P)H] small sub.	0.81	0.2	0.2
	CDW43_RS11595	Nitrite reductase [NAD(P)H] large sub.	0.84	0.4	0.3
	Other genes involved in nitrogen pathway
	CDW43_RS01710	fnr transcriptional regulator	1.41	3.9	5.5
	CDW43_RS00310	nnrS protein involved in response to NO	2.71	4.4	12.2
	CDW43_RS01960	nnrS	0.99	7.4	7.3
	CDW43_RS08605	nnrS	1.61	3.8	6.4
	CDW43_RS00315	nsrR transcriptional regulator	1.90	17.4	33.8
	CDW43_RS01825	nsrR	0.86	2.3	2.0
	CDW43_RS00320	dnrN /norA/yftE NO-dependent regulator	1.98	25.0	48.7
	CDW43_RS01015	Ammonium transporter	1.01	1.0	1.0
	CDW43_RS05605	Ammonium transporter	0.48	0.9	0.5
	CDW43_RS01835	Nitric oxide dioxygenase	0.79	2.7	2.2
	CDW43_RS09635	Nitric oxide dioxygenase	0.65	8.0	5.0
	CDW43_RS01010	Nitrogen regulatory protein P-II	0.85	1.2	1.0
	CDW43_RS06615	Nitrogen regulatory protein P-II	0.80	7.0	5.5
	CDW43_RS13530	Nitrogen regulatory protein P-II	0.74	1.2	0.8
	CDW43_RS01095	Nitrogen regulation protein NtrY	1.30	1.4	1.9
	CDW43_RS01100	Nitrogen regulation protein NtrX	1.35	3.6	4.9
	CDW43_RS06315	Nitrogen regulation protein NtrB	0.60	1.1	0.6
				(continued on next page)

Table 4

 4 (continued) Transcript levels relative to the expression of the methanol dehydrogenase small unit (mxaI ) set to 100 (see Data S6). For polycistronique operons, the values are the average values of individual genes. na, not applicable. Locus tag refers to the genome sequence of strain GP59 (GenBank accession number CP021973). Gene/function is based on RAST annotations, see Data S4.

	Locus Tag	Gene/function	a Expression ratio	b mxaI = 100
			JAM1/GP59	GP59	JAM1
	CDW43_RS06320	Nitrogen regulation protein NR(I)	1.12	0.7	0.8
	CDW43_RS15625	PTS IIA-like nitrogen-regulator PtsN	0.95	12.1	11.1
	Notes.
	a Ratio of JAM1 transcripts per million (TPM) by GP59 TPM.

b

Table 5 Relative expression of selected genes under anoxic culture conditions.

 5 Gene expression is relative to mxaI set to 100 (2602 TPM for strain GP59 and 4060 TPM for strain JAM1). Locus tag: from GenBank accession number CP021973. Functions are based on RAST annotations.

	Locus tag	Function

Geoffroy et al. (2018), PeerJ, DOI 10.7717/peerj.4679

Geoffroy et al. (2018), PeerJ, DOI 10.7717/peerj.4679 9/25

Geoffroy et al. (2018), PeerJ, DOI 10.7717/peerj.4679 20/25

Geoffroy et al. (2018), PeerJ, DOI 10.7717/peerj.4679 25/25

ACKNOWLEDGEMENTS

We thank Karla Vasquez for her technical assistance.

Data Availability

The following information was supplied regarding data availability:

The raw data and supplemental documents are available as Supplemental Files.

ADDITIONAL INFORMATION AND DECLARATIONS Funding

This research was supported by a grant to Richard Villemur from the Natural Sciences and Engineering Research Council of Canada # RGPIN-2016-06061. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures

The following grant information was disclosed by the authors: Natural Sciences and Engineering Research Council of Canada: #RGPIN-2016-06061.

Competing Interests

The authors declare there are no competing interests.

Author Contributions

• Valérie Geoffroy and Geneviève Payette conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, approved the final draft.

• Florian Mauffrey conceived and designed the experiments, performed the experiments, analyzed the data, authored or reviewed drafts of the paper.

Supplemental Information

Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.4679#supplemental-information.


Etta Grover-Silva

Miguel Heleno

Salman Mashayekh

Gonçalo Cardoso

Robin Girard

George Kariniotakis

A Stochastic Optimal Power Flow for Scheduling Flexible Resources in Microgrids Operation

Keywords: demand response, microgrids, optimal power flow, photovoltaics, stochastic optimization, storage

published or not. The documents may come   L'archive ouverte pluridisciplinaire

A Stochastic Optimal Power Flow for Scheduling Flexible Resources in Microgrids Operation

Etta Grover-Silva1,3 , Miguel Heleno2 , Salman Mashayekh 2 , Gonçalo Cardoso 2 , Robin Girard

I. INTRODUCTION

t the distribution grid level, uncertainties in renewable generation and load consumption represent a challenge to network operation, namely for day ahead planning of Distributed Energy Resources (DERs), such as grid connected storage, controllable loads or photovoltaic (PV) control strategies, implemented in real time by a distribution management system (DMS). These challenges are magnified in microgrids, where uncertainties are higher due to minimal aggregation and smoothing effects. Since microgrids are more easily perturbed by DERs, an accurate control is needed to manage multiple electric storage systems, load devices and generation units, while ensuring a stable and reliable operation of the microgrid network and minimizing costs [START_REF] Aris | Microgrids Control Issues in Microgrids: Architectures and Control[END_REF] [START_REF] Mashayekh | A mixed integer linear programming approach for optimal DER portfolio, sizing, and placement in multi-energy microgrids[END_REF].

Due to high uncertainties in load and renewable generation, microgrid control requires advanced forecasting tools and robust scheduling of controllable devices to guarantee power quality and security of supply. In particular, the control of individual loads, e.g. heating, ventilation and air-conditioning (HVAC) systems [START_REF] Alqurashi | Treatment of uncertainty for next generation power systems: State-of-the-art in stochastic optimization[END_REF], brings new sources of uncertainty to the day ahead planning of DERs, such as ambient temperature, building occupancy and consumption habits. This uncertainty has a modest impact on grid operations when aggregated at the distribution level but it becomes relevant at the microgrid scale where a finer control is needed.

Optimization algorithms have been presented in the literature to solve the problem of day ahead scheduling of microgrid dispatchable resources. Numerous examples of deterministic [START_REF] Maulik | Optimal operation of microgrid using four different optimization techniques[END_REF] [START_REF] Baldi | Joint energy demand and thermal comfort optimization in photovoltaic-equipped interconnected microgrids[END_REF], stochastic [START_REF] Wang | Coordinated Energy Management of Networked Microgrids in Distribution Systems[END_REF] [START_REF] Good | Optimization Under Uncertainty of Thermal Storage-Based Flexible Demand Response With Quantification of Residential Users" Discomfort[END_REF] and hybrid [START_REF] Wang | Integrated scheduling of energy supply and demand in microgrids under uncertainty: A robust multi-objective optimization approach[END_REF] [START_REF] Krishnanand | Optimal appliance scheduling in building operating systems for cost-effective energy management[END_REF] approaches to plan and operate renewable intensive microgrids are presented in the literature. Optimization methods include quadratic programming (QP) [START_REF] Grover-Silva | Stochastic Multi-temporal Optimal Power Flow Approach for the Management of Grid Connected Storage[END_REF], as well as heuristic and meta-heuristic techniques [START_REF] Maulik | Optimal operation of microgrid using four different optimization techniques[END_REF] [START_REF] Gazijahani | Optimal day ahead power scheduling of microgrids considering demand and generation uncertainties[END_REF]. In the optimal scheduling of DERs in multi-node microgrids, heuristics have the advantage of enabling exact network constraints [START_REF] Maulik | Optimal operation of microgrid using four different optimization techniques[END_REF], while QP requires a convex relaxation of power flow equations [START_REF] Krishnanand | Optimal appliance scheduling in building operating systems for cost-effective energy management[END_REF]. Due to the random aspect of search techniques in heuristic methods, calculation time can be high and the global optimal is not guaranteed. However, QP methods perform significantly better in terms of computational time. Thus, when combined with techniques that guarantee accuracy of the power flow calculations, e. g. linear cuts [START_REF] Abdelouadoud | Optimal power flow of a distribution system based on increasingly tight cutting planes added to a second order cone relaxation[END_REF], they become a better solution.

Stochastic approaches have been used in optimal operation of microgrids to capture uncertainties of renewable sources [START_REF] Su | Stochastic Energy Scheduling in Microgrids With Intermittent Renewable Energy Resources[END_REF]. Primarily, these strategies include either scenario trees [START_REF] Trovato | Role and Benefits of Flexible Thermostatically Controlled Loads in Future Low-Carbon Systems[END_REF] [START_REF] Good | Optimization Under Uncertainty of Thermal Storage-Based Flexible Demand Response With Quantification of Residential Users" Discomfort[END_REF] or statistical parameters of the stochastic variables [START_REF] Nikmehr | Probabilistic optimal scheduling of networked microgrids considering time-based demand response programs under uncertainty[END_REF] [11] [START_REF] Zhang | Robust Energy Management for Microgrids With High-Penetration Renewables[END_REF] that are integrated into the optimization problem.

A

Monte Carlo simulation along with the distribution functions for generators and load are used to generate scenarios in [START_REF] Nikmehr | Probabilistic optimal scheduling of networked microgrids considering time-based demand response programs under uncertainty[END_REF]. Scenarios are constructed by analyzing the mean, standard deviation and probability density functions of load and generation in [START_REF] Gazijahani | Optimal day ahead power scheduling of microgrids considering demand and generation uncertainties[END_REF]. Upper and lower bounds on generation and load are considered in [START_REF] Zhang | Robust Energy Management for Microgrids With High-Penetration Renewables[END_REF]. A scenario tree is developed to represent stochastic variables such as temperature, electricity prices and consumer occupancy through the calculation of quantiles and consideration of the probability density function (PDF) of historical data [START_REF] Good | Optimization Under Uncertainty of Thermal Storage-Based Flexible Demand Response With Quantification of Residential Users" Discomfort[END_REF].

The day ahead operation of microgrids includes optimal scheduling of multiple DER technologies. Besides the generation and storage control solutions, demand response (DR) has been a valuable resource to compensate the variability of the renewable sources, especially through the control of thermal loads, such as HVAC and Electric Water Heaters (EWH). In fact, as shown in [START_REF] Wang | Integrated scheduling of energy supply and demand in microgrids under uncertainty: A robust multi-objective optimization approach[END_REF], load control can significantly reduce microgrid operation costs as well as CO 2 emissions. Two primary modeling strategies are presented in the literature for DR consideration: an aggregated model or individual modeling of devices. Aggregated models make acceptable assumptions about individual devices [START_REF] Trovato | Role and Benefits of Flexible Thermostatically Controlled Loads in Future Low-Carbon Systems[END_REF] and improve aggregated controllability of the microgrid, but the comfort of individual end users is not modeled in detail. Thus, individual load models become more appropriate for small scale applications (e.g. buildings) where a detailed comfort representation is required. In [START_REF] Krishnanand | Optimal appliance scheduling in building operating systems for cost-effective energy management[END_REF], individual load models are used in optimization of building operations with DR. A deterministic approach that considers end-user comfort constraints and PV for a 3 building micro-grid is detailed in [START_REF] Baldi | Joint energy demand and thermal comfort optimization in photovoltaic-equipped interconnected microgrids[END_REF]. An algorithm proposing an economic penalty for violations in thermal comfort constraints is presented in [START_REF] Good | Optimization Under Uncertainty of Thermal Storage-Based Flexible Demand Response With Quantification of Residential Users" Discomfort[END_REF] however, this paper does not consider the electric network and instead performs only an energy balance.

A majority of the mentioned citations take into account the losses in the electrical lines in a two-step process and do not integrate a full AC optimal power flow (AC-OPF) into the optimization problem [START_REF] Su | Stochastic Energy Scheduling in Microgrids With Intermittent Renewable Energy Resources[END_REF] [11] [START_REF] Nikmehr | Probabilistic optimal scheduling of networked microgrids considering time-based demand response programs under uncertainty[END_REF]. This paper presents a novel method for day ahead scheduling of loads and DERs that has a low calculation burden while considering network constraints. To the authors knowledge, it is the first time that a full AC-OPF algorithm is used while considering thermal comfort constraints of end users. Moreover, the presented model adds on recent innovations in the field of stochastic AC-OPF [START_REF] Baradar | A stochastic SOCP optimal power flow with wind power uncertainty[END_REF], by expanding the formulation to accommodate multiple sources of uncertainty relevant for microgrid operations: PV generation, load and ambient temperature.

The algorithm uses the second order cone program (SOCP) convex relaxation of the power flow equations proposed in [START_REF] Abdelouadoud | Optimal power flow of a distribution system based on increasingly tight cutting planes added to a second order cone relaxation[END_REF] to ensure accuracy in load flow constraints with a relatively low computational burden. Also, to avoid relaxations to be inexact, especially in periods of high DER injection, linear cuts are added to the problem to guarantee exactness on load flow constraints.

The microgrid flexible resources considered in the optimal scheduling model are electrochemical storage and thermal loads, namely space heating/cooling and EWH devices. To include thermal loads in the formulation, this paper integrates into the stochastic OPF different developments on individual load models [START_REF] Krishnanand | Optimal appliance scheduling in building operating systems for cost-effective energy management[END_REF], thermal comfort constraints [START_REF] Heleno | Availability and flexibility of loads for the provision of reserve[END_REF] and uncertainty in ambient temperature [START_REF] Good | Optimization Under Uncertainty of Thermal Storage-Based Flexible Demand Response With Quantification of Residential Users" Discomfort[END_REF]. The result is a suitable DER scheduling method, based on one of the most detailed formulation of its kindcombining first-order thermal models with second order security constraintsto support control at the microgrid level. Thus, the main contributions of the paper are the following:  a multi-period SOCP adapted to consider uncertainties through scenarios of generation, load, hot water consumption and ambient temperature to account for thermal loads;  the optimal day-ahead scheduling of microgrid flexibilities, considering grid constraints, end-user comfort constraints, and the multi-temporal dispatch of different DER technologies;  the behind-the-meter individual loads devices modeling and scheduling for optimal DR strategies, constrained by the comfort of end-users, and integrated with the microgrid stochastic dispatch. By combining these contributions in a single stochastic AC-OPF, the authors aim at providing a valuable discussion on the implications of generation and load uncertainties for microgrids control and the resulting effects on end-user comfort while considering demand side management. Following this introductory section, section II describes the novel formulation introduced by this paper. Section III describes a case study to demonstrate the utility of the day-ahead scheduling strategies produced, and section IV discusses final conclusions.

II. STOCHASTIC OPTIMAL POWER FLOW METHOD

A. Methodology Overview

This section proposes a multi period stochastic optimization method for the day ahead scheduling of DERs in a microgrid. The DERs considered are electric storage and controllable (CL) buildings" water and space heating/cooling loads, such as EWH and commercial HVAC units. This stochastic approach considers uncertainty in the baseline uncontrollable loads (UL) (such as lighting, cooking appliances, electronic devices and phantom loads), PV generation, ambient temperature and hot water consumption. These uncertainties are considered in the form of forecast scenarios which are generated from probabilistic forecasts taking into account the spatial and temporal correlations in the processes. High and low scenarios for each variable are selected through a scenario reduction strategy and are assumed equally probable as described later in section II C.

The benefit introduced by the stochastic approach is measured by the value of the stochastic solution (VSS).This consists of comparing the expected value of perfect information (EVPI) given the stochastic solution and the deterministic average solution [START_REF] Eajal | Stochastic Centralized Dispatch Scheme for AC/DC Hybrid Smart Distribution Systems[END_REF].A schematic showing the methodology as a flow chart is shown in Fig. 1.

B. Formulation

The proposed formulation (1)-( 22) is a multi-period (t) multi-scenario (s) optimal power flow that aims at reducing the day ahead microgrid operation costs through scheduling of batteries (𝑃 st ,j,t ) and controllable thermal loads, EWH (𝑃 ewh ,j,t ) and HVAC (𝑃 hvac ,j,t ), located at the nodes (j) of the microgrid.

1) Objective function

The operation cost function (1) considers differentiated rates for energy imports and exports, following the current regulatory mechanisms adopted by several European countries to promote self-consumption. Hence, the energy exports at the point of common coupling are remunerated at wholesale electricity market price, while the energy consumption costs are charged at the final electricity price, which corresponds to the hourly electricity market price with fixed rates, e. g. due to transmission and distribution cost. The comfort constraints of space heating systems are also considered in the objective (1)

2) Power flow constraints

The constraints of the problem include the nodal power balance considering different DER units (2)-( 5). Equation ( 6) describes the convex relaxations of the line constraints. The result of each OPF calculation for each time step is compared to a forward backward sweep power flow calculation to verify that the convex relaxation of the line constraint equation ( 6) is exact. If the solution is not exact, linear cuts are added to the problem to guarantee exactness as explained in [START_REF] Abdelouadoud | Optimal power flow of a distribution system based on increasingly tight cutting planes added to a second order cone relaxation[END_REF].

𝑃 ij
𝜗 j,s,t = 𝜗 i,s,t -2 r ij 𝑃 ij,s,t + x ij 𝑄 ij,s,t + r ij 2 + x ij 2 ℓ ij,s,t (3)
V 2 ≤ 𝜗 j,s,t ≤ V 2 (5)

ℓ ij,s,t ≥ 𝑃 ij,s,t 2 + 𝑄 ij,s,t 2 𝜗 i,s,t (6)
3) Battery system constraints Equations ( 7)- [START_REF] Grover-Silva | Stochastic Multi-temporal Optimal Power Flow Approach for the Management of Grid Connected Storage[END_REF] represent the battery limits regarding power and state of charge.

-P st ,j ≤ 𝑃 st ,j,t ≤ P st ,j

soc st ,j ≤ 𝑠𝑜𝑐 st ,j,t ≤ soc st,j

𝑃 st ,j,t = 𝑃 st +,j,t + 𝑃 st -,j,t

𝑠𝑜𝑐 st ,j,t = 𝑠𝑜𝑐 st ,j,t-1 + tη st 𝑃 st +,j,t + t η st 𝑃 st -,j,t (10)

4) Thermal comfort constraints

The thermal comfort constraints associated with the individual HVAC and EWH controllable devices are shown in ( 12)- [START_REF] Eajal | Stochastic Centralized Dispatch Scheme for AC/DC Hybrid Smart Distribution Systems[END_REF]. The division of over and under heating in equation ( 20)-( 22) allows for a piecewise linear penalty function of thermal constraint violations in the objective function. 𝑃 cl,j,t = 𝑃 hvac ,j,t + 𝑃 ewh ,j,t

𝑃 ewh ,d,t 𝑑 = 𝑃 ewh ,j,t

0 ≤ 𝑃 ewh ,d,t ≤ P ewh ,d (13)
θ hvac ,d -𝜃 hvac ,d,t = ∆𝜃 low ,d,t (20)
𝜃 hvac ,d,t -θ hvac ,d = ∆𝜃 high ,d,t (21)
These thermal equations are the first order physically-based load modesconsidering the thermal capacity (C), resistance (R), and heat loss constant (α) -to describe the temperature behavior of thermal systems.

5) Stochastic variable and scenario generation

Controllable variables include the active power of EWH for residential buildings and HVAC thermal loads for commercial buildings (𝑃 ewh ,d,t and 𝑃 hvac ,d,t ) and the active power of battery systems (𝑃 st ,j,t ).On the other hand, the baseline uncontrollable load (𝑃 ul ,j,t ) is also is considered as a state variable. A table summarizing the controllable variables, stochastic variables and scenario dependent variables is found in Table I.

The stochastic variables are represented in the linear constraint matrix of the optimization problem through parallel multi-period scenarios. This technique to integrate uncertainties into and OPF is classified as a probabilistic scenario-based technique for taking into account uncertainties in power systems as classified in [START_REF] Aien | A somprehensive review on uncertainty modeling techniques in power system studies[END_REF] a.k.a. a deterministic equivalent formulation of the stochastic problem. Here we apply the approach proposed in [START_REF] Bakirtzis | Multiple time resolution stochastic scheduling for systems with high renewable penetration[END_REF], where a similar scenario representation in an OPF scheduling algorithm has been implemented. Two main steps are necessary, the generation of scenarios and the reduction of scenarios. That paper uses a Monte Carlo scenario generation method applied to wind turbine generation and load profiles based on [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF]. The scenario reduction technique that is used is based on probabilistic distance and fast-forward selection as described in [START_REF] Conejo | Decision Making Under Uncertainty in Electricity Markets[END_REF].

Here we use the same scenario generation technique as [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF] for the UL, ambient temperature and EWH scenarios. A three-month historical period is used to calculate the quantiles and the covariance matrix to generate normal Gaussian scenarios.

This technique is less effective when applied to PV production scenarios due to the fact that PV production profile has a very strong correlation associated with the irradiation which depends on the course of the sun during the day. This strong correlation with irradiation may dilute the other causes of variation in PV production such as cloud cover. The PV production profiles were therefore normalized by the clear sky index before applying the scenario generation method. This allows for a more precise analysis of inter-temporal variation due to cloud cover or other phenomena that are not correlated with irradiation. For further discussion on the necessity to stationnarise PV production time-series when modeling spatio-temporal correlations and more sophisticated stationnarisation techniques for that purpose, refer to [START_REF] Agoua | Short-Term Spatio-Temporal Forecasting of Photovoltaic Power Production[END_REF].

It is noted that a limitation of the approach in [START_REF] Bakirtzis | Multiple time resolution stochastic scheduling for systems with high renewable penetration[END_REF] concerns the fact that the scenario generation techniques does not consider possible correlations between the stochastic variables considered. This could be part of the perspectives of the current work. The stochastic OPF optimization, through the use of multiple scenarios as input, aims at finding a solution that provides "hedging" to the considered physical system against uncertainties. This means that the system is prepared to face more situations than when optimized through the simpler deterministic approach. This strategy may consequently involve higher costs. To be able to evaluate this risk for higher costs we have adopted a simple but intuitive scenario reduction method that is based on choosing the "extreme" upper and lower scenarios resulting from the scenario generation step and then combine these opposite situations that these scenarios reflect. This can be considered as a pessimistic approach that could lead to amplified "hedging" costs compared to the deterministic case. Depending on the system this does not necessarily mean that there is no margin to reduce costs through the deterministic approach. It is a matter of tradeoff between the hedging cost mentioned above and the impact of the deviations from a deterministic schedule. The overall approach remains though generic as one can replace by scenarios resulting from more sophisticated reduction method (i.e. from in [START_REF] Conejo | Decision Making Under Uncertainty in Electricity Markets[END_REF]) as done by the authors in [START_REF] Bakirtzis | Multiple time resolution stochastic scheduling for systems with high renewable penetration[END_REF]. More precisely, here, three scenarios were selected based on the total cumulative values of the day. The three scenarios are chosen by selecting the maximum, minimum and closest to average value of the cumulative values in order to produce a high, low and average scenario for each variable.

C. Performance Evaluation of Stochastic Method

The stochastic approach is evaluated through VSS and EVPI. To calculate these performance indicators, the day-ahead schedules obtained in the deterministic approach and those obtained in the stochastic approach are benchmarked against the actual observation of day-of conditions.VSS is the difference between the stochastic method and the deterministic one. The expected value of the actual observation, EVPI, is the absolute difference between the expected value with the data of the actual observation, and the expected value without the actual observationeither the stochastic solution or the deterministic solution. When evaluating each of these cases, the set points for EWH, HVAC and battery power are implemented with no intra-day adjustments for the real conditions. The thermal equations are used to simulate the evolution of temperatures in the buildings and in the hot water tanks. A forward backward sweep power flow calculation is performed to calculate the current and voltage of each node at each time step. Energy costs, grid constraint violations and thermal comfort profiles are then analyzed to assess the comfort of end-users in comparison with the economic performance of the optimization strategies.

III. CASE STUDY

In this section, the proposed methodology presented above is applied to a realistic community microgrid, located in France in the area of Grenoble city. To facilitate the interpretation and to allow the replicability of the results, the microgrid network is assumed to be the well-known IEEE 37 nodes distribution circuit(Fig. 2). However, generation and load were modified as described below.

A. Generation and Load Data

Load profiles including EWH, HVAC and uncontrollable loads are generated using a bottom up load simulator detailed in [START_REF] Barbier | A novel approach for electric load curve holistic modeling and simulation[END_REF]. This simulator produces a group of individual commercial and residential building load profiles. These profiles are generated to be statistically accurate representations of residential and commercial customer proportion, electric heating, building surface area and population using the INSEE building inventory database of France, and distributed randomly across the network. The medium voltage feeder is assumed to have a 5 MVA transformer serving 5 low voltage substations for a total of 312 clients, of which 300 are residential and 12 are commercial. Of the 300 residential clients, 49 residential hot water heaters are controllable. It was assumed that all 12 commercial clients have controllable HVAC systems. In addition, a total of 1.2 MW uncontrollable load, 155 kW of controllable EWH and 308 kW of controllable HVAC are considered.

A total capacity of 2.62 MW of PV is distributed over 5 nodes, with production curves based on a real PV plant in Grenoble [START_REF] Solaire Lgm | Solaire lgm[END_REF], normalized by the nominal power installed in each node. In addition, all 5 PV nodes are assumed to have a battery system with 250 kW nominal power and 500 kWh nominal capacity, totaling 1.25 MW and 2.5 MWh. The UL, CL and DER characteristics for each node can be found in Table II. The parameters used for the HVAC and EWH units are as follows: for EWH, the maximum power per device is between 2.0-6.0 kW, thermal capacities and heat loss coefficients are within 0.

B. Results

The results presented below focus on three primary topics: i) economic benefit of using a stochastic approach over a deterministic one, ii) economic benefits of combining stochastically managed storage devices with controllable loads and iii) thermal comfort improvements with stochastic techniques.

The annual operational costs of the microgrid are calculated using both the stochastic and deterministic day-ahead scheduling strategies. The deterministic case (denoted as "Det" in the figures) uses the average daily cumulative value of the available forecasts for all stochastic variables. The stochastic case (denoted as "St" in the figures) uses combinations of high and low forecasts for each stochastic variable shown in Table I.

Multiple case studies are tested to quantify the effect of each stochastic variable on the total annual operational cost. The cases consider high (H) and low (L) scenarios and combine different stochastic variables to analyze the effects of each stochastic variable independently as well as their compounded effects. Cases when a stochastic variable is not being evaluated the average scenario is used (M). Table III details the labels for each case study evaluated in this section. Fig. 3 shows the stochastic performance measures of the annual operational cost. The case study with CL shows that stochastic strategies result in lower annual costs. However, no major economic advantage is seen in using 4 scenarios instead of 2. In fact, the 2 initial scenarios already represent extreme conditions of PV and UL and they have a dominant impact in the stochastic optimization.

Fig. 3. Annual VSS (top) and EVPI (bottom) costs of case studies

As seen in Fig. 3, when integrating controllable load into the optimal scheduling problem, annual costs are reduced. This reduction is primarily due to shifting HVAC and EWH to less expensive periods.

Figure 4 shows the daily operation of the microgrid, presenting a comparison of costs, battery and load scheduling between deterministic and stochastic approaches. The stochastic day ahead strategy significantly improves the operation costs of the microgrid system when simulated with real day of conditions. This improvement can be explained by the observation of the microgrid operation during the transition between hours 6 to 9 and 14 to 19. In these periods deterministic approach is unable to respond accurately to situations where PV production deviates from the predicted average. For example, the deterministic scheduling solution suggests discharging the battery between 6 and 9 AM, preparing it to store the PV surplus during the expected sunny day. However, the observed PV generation is significantly lower than expected and this decision has a negative economic impact in the subsequent periods. In contrast, since the stochastic approach considers from the beginning the scenario of low PV generation, it suggests a more conservative discharge during the early morning. Instead, the battery is significantly discharged in the period10-13 hours, where the electricity price is higher.

From the analysis performed, a tradeoff was observed between annual operational costs and thermal comfort of the users. This tradeoff can be represented as a Pareto-optimality state with annual operational costs and end user comfort defining the Pareto-frontier.

To evaluate the trade offs of annual operational costs and user comfort additional scenarios are considered. These scenarios include high and low ambiant temperature scenarios and high and low hot water consumption scenarios. The cases "St w T 1 €" and "St w T 10 €" associates a 1 €/°Ch or 10 €/°Ch penalty for comfort constraint violations. The high and low temperature scenarios with economic penalties results in a more robust management of the heating and hot water loads. This robust management tends to keep the temperatures in a middle range as opposed to only heating the minimal amount. Therefore, increasing annual operational costs and improving user comfort. As shown in Fig. 5, considering uncertainties in ambient temperature is a Pareto improvement for end user comfort in HVAC loads while not considering these uncertainties is a Pareto improvement for annual costs. The value chosen for the cost penalties, 1 €/°Ch and 10 €/°Ch, effected the optimization differently. With higher cost penalties for temperature violations, annual costs were higher but fewer comfort violations resulted. In all cases the stochastic algorithm results in lower annual operational costs for the same number of comfort constraint violations when comparing stochastic and deterministic approaches.

C. Algorithm Performance

The algorithms presented in this paper have been implemented in Python and solved using the MOSEK SOCP solver on an 8-core, 3.4 GHz CPU. Due to the fact that the stochastic algorithm takes into account multiple scenarios the calculation time of this algorithm is higher. A performance analysis was completed to compare the time of calculation for the deterministic algorithm and varying amounts of scenarios in the stochastic algorithm. The average calculation time of stochastic and deterministic methods for 24 coupled time steps is shown in Table IV. Therefore, the stochastic algorithm is about 3 times slower with the consideration of 2 scenarios and 6 times slower with the consideration of 4 scenarios.

IV. CONCLUSION

This paper presents a multi-temporal stochastic method that performs a centralized day-ahead economic scheduling of storage and controllable loads in a microgrid. The method considers network constraints using a SOCP approach in order to guarantee the security of the microgrid operation. Moreover, end-user comfort is ensured during the optimal scheduling of thermal loads, by considering first order thermal constraints. Finally, the stochastic formulation considers different scenarios, generated based on a Monte Carlo technique that takes into account spatio-temporal correlations of the variables.

The use of a stochastic approach resulted in a reduction of the microgrid operation costs in comparison with the deterministic strategy, especially in periods where weather conditions and baseline load deviate from the average. The results also show that significant savings can be achieved by harnessing thermal controllable loads in day-ahead scheduling of microgrids. In particular, considering detailed physical models and end user comfort constraints is important to effectively implement load scheduling without affecting end user comfort. However, especially in stochastic applications, meeting comfort constraints in all possible scenarios can lead to extreme conservative scheduling solutions with higher costs. Therefore, formulations based on comfort penalties, also explored in paper, prove to be suitable alternatives for this kind of constraint. Smartgrids Group. He has authored more than 220 scientific publications in journals and conferences. He has been involved as participant or coordinator in more than 40 R&D projects in the fields of renewable energies and distributed generation.

Fig. 1 .

 1 Fig. 1. Flow chart of proposed methodology

𝑚𝑖𝑛

  0877 -0.2925 kWh/ o C and 0.0004-0.0012 kW/ o C, respectively. Cold water intake temperature (θ int ) and usage temperatures (θ out ) are12 o C and65 o C, while temperature limits are between 60 o C and 80 o C. For individual buildings HVAC systems, the maximum power is between 2.44 -158.67, C and R values are within 0.2244 -1318.4959 kWh/ o C and 0.0127-21.0012 o C/kW, respectively. The comfort temperatures are between19 o C and 26 o C.The cost of discomfort for under heating and overheating was considered to be

Fig. 2 .

 2 Fig. 2. Medium voltage microgrid case study feeder. C indicates a node with customer load, PV + B indicates nodes with PV and batteries

Fig. 4 .

 4 Fig. 4. Example day where stochastic scheduling results in lower operational costs than the deterministic one. From top to bottom: cost of electricity (1), PV scenarios and real PV production (2), uncontrollable load (UL) scenarios and real UL (3), stochastic and deterministic battery schedule (4), stochastic and deterministic controllable load schedule (5).

Fig. 5 .

 5 Fig. 5. EVPI vs comfort constraint violations of HVAC and EWH systems

  . c c,t 𝑡𝑃 0+,s,t + c e,t 𝑡𝑃 0-,s,t + c cf ∆𝜃 low ,d,t + ∆𝜃 high ,d,t

	s	t	d

  𝜃 ewh ,d,t-1 -𝜃 int ,d,t -ν d,t c w θ out -θ in + 𝑃 ewh ,d,t

	𝜃 ewh ,d,t = 𝜃 ewh ,d,t-1 +	t C d	-α d (16)
	𝑑	𝑃 hvac ,d,t	= 𝑃 hvac ,j,t		(17)
	0 ≤ 𝑃 hvac ,d,t ≤ P hvac ,d		(18)
	𝜃 hvac ,d,t = 𝜃 hvac ,d,t-1 -	t C d R d	𝜃 hvac ,d,t-1 -θ ext ,t + η d R d 𝑃 hvac ,d,t	(19)
	θ hvac ,d -∆θ low ,d ≤ 𝜃 hvac ,d,t ≤ θ hvac ,d + ∆θ high ,d
	θ w < 𝜃 ewh ,d,t < θ w		(15)

Table . 1

 . Variable types TypeApplicable Variables Controllable Variables P st (P st+ , P st-), P cl (P hvac ,

		P ewh )
	Scenario Dependent	P ij , ℓ ij , Q ij , ϑ i
	Variables
	Stochastic Variables	P pv , P ul , v d,t , θ ext

Table . 3

 . Case

			study labels
	Scenario	ST NO CL 2 S	ST NO CL 4 S	ST W T 10 €	ST W T 1 €
	1 2 3	P 𝑝𝑣 𝐿 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐿 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 𝑣 -	P 𝑝𝑣 𝐿 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐿 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀	P 𝑝𝑣 𝐿 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝐻 θ 𝑒𝑥𝑡 𝐿 P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐿 P 𝑒𝑤 ℎ 𝐿 θ 𝑒𝑥𝑡 𝐻 P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝐻 θ 𝑒𝑥𝑡 𝐻	P 𝑝𝑣 𝐿 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝐻 θ 𝑒𝑥𝑡 𝐿 P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐿 θ 𝑒𝑥𝑡 𝐻 𝐿 P 𝑒𝑤 ℎ P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝐻 θ 𝑒𝑥𝑡 𝐻
	4	-	P pv L P ul L P ewh M θ ext M	P pv L P ul L P ewh L θ ext L	P pv L P ul L P ewh L θ ext L
	Scenario	ST 2 S	ST C	ST PV	ST 4 S
	1 2 3	P 𝑝𝑣 𝐿 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐿 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 -	P 𝑝𝑣 𝑀 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 P 𝑝𝑣 𝑀 P 𝑢𝑙 𝐿 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 -	P 𝑝𝑣 𝐻 P 𝑢𝑙 𝑀 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 P 𝑝𝑣 𝐿 P 𝑢𝑙 𝑀 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 -	P 𝑝𝑣 𝐿 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀 P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐿 P 𝑒𝑤 ℎ 𝑀 𝑀 θ 𝑒𝑥𝑡 P 𝑝𝑣 𝐻 P 𝑢𝑙 𝐻 P 𝑒𝑤 ℎ 𝑀 θ 𝑒𝑥𝑡 𝑀
	4	-	-	-	P pv L P ul L P ewh M θ ext M

Table . 4

 . Algorithmic performance

		Number	Coupled	CL	Time
	Analysis	of	Timesteps	present	(s)
		Scenarios
	Deterministic	1	24	No	3
	Stochastic	2	24	No	8
	Stochastic	4	24	No	14
	Deterministic	1	24	Yes	7
	Stochastic	2	24	Yes	15
	Stochastic	4	24	Yes	45

MINES ParisTech, PSL Research University, PERSEE -Center for Processes, Renewable Energies and Energy Systems (e-mail: etta.groversilva@gmail.com, robin.girard@mines-paristech.fr, georges.kariniotakis@mines-paristech.fr)

Emest Orlando Lawrence Berkeley National Laboratory (LBNL) (email: miguelheleno, smashayekh, gfcardoso @lbl.gov)

ADEME -French Environment and Energy Management Agency -France

V. ACKNOWLEDGEMENTS The work is performed in the frame of the project SENSIBLE (Grant No 645963) funded in part by the European Commission under the Horizon 2020 Framework Program as well as a PhD scholarship grant supported by the French Energy Agency ADEME. The international collaboration between MINES ParisTech PERSEE and LBNL was made possible by a mobility grant in the frame of the Electra IRP project (Grant No 609687) funded in part by the European Commission under the 7 th Framework Program and of the MINES ParisTech Foundation. This work is also partially funded by the Office of Electricity Delivery and Energy Reliability, Smart Grid R&D Program of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.