DNA Methylation Epigenetically Regulates Gene Expression in Burkholderia cenocepacia and Controls Biofilm Formation, Cell Aggregation, and Motility

CF patients diagnosed with Burkholderia cenocepacia infections often experience rapid deterioration of lung function, known as cepacia syndrome. B. cenocepacia has a large multireplicon genome, and much remains to be learned about regulation of gene expression in this organism. From studies in other (model) organisms, it is known that epigenetic changes through DNA methylation play an important role in this regulation. The identification of B. cenocepacia genes of which the expression is regulated by DNA methylation and identification of the regulatory systems involved in this methylation are likely to advance the biological understanding of B. cenocepacia cell adaptation via epigenetic regulation. In time, this might lead to novel approaches to tackle B. cenocepacia infections in CF patients.

accompanied by bacteremia and sepsis. If left untreated, cepacia syndrome can lead to death within weeks (8,9). The genome of B. cenocepacia is complex (with usually three large replicons), with a high GC content (67%) and large size, comprising approximately 8.06 Mb (10). The species has been classified into different phylogenetic clusters and subdivided into lineages, including the highly transmissible ET-12 lineage that harbors B. cenocepacia strains J2315 and K56-2 (11,12).
Epigenetics is the study of heritable changes in gene expression without changes in the actual genomic sequence. In bacterial genomes, epigenetic control is exerted by DNA methyltransferase enzymes (MTases) (13)(14)(15). DNA MTases originate from restriction-modification (RM) systems, early defense mechanisms in bacteria with an active interplay between endonucleases and DNA MTases, which cleave foreign DNA but protect the organism's own genome. In addition, discovery of orphan DNA MTases, enzymes without a cognate endonuclease, shows that DNA MTases are not exclusively dependent on the presence of the restriction part to function as regulator of gene expression (16).
DNA MTases interact with specific DNA recognition sites and transfer a CH 3 group from a methyl donor, mostly S-adenosylmethionine (SAM), to a cytosine (C 5 -methyl cytosine or N 4 -methyl cytosine) or adenine (N 6 -methyl adenine) base (17,18). As methylated bases change the binding affinity of DNA binding proteins, methylation at regulatory regions allows bacteria to regulate gene expression at the level of transcription (19,20). While both cytosine and adenine methylation occur in eukaryotic and prokaryotic cells, C 5 -methyl cytosine is the archetypal eukaryotic base methylation signature (16,21). Conversely, in prokaryotes, N 6 -methyl adenine is the most important base modification involved in gene expression regulation (22). In addition to this, studies with DNA MTases Dam (deoxyadenosine methyltransferase) and Dcm (DNA cytosine MTase) in Escherichia coli have demonstrated that, besides having a regulatory function, DNA MTases also take part in crucial cellular processes like DNA replication initiation or methyl-directed mismatch repair (21,23).
Detection of (genome-wide) DNA methylation patterns has been challenging in the past. The use of specific restriction enzymes with affinity for methylated sites, followed by a comparison of the resulting fragment lengths, gave a good impression of methylation of the treated DNA at one particular area, but global methylation analysis was until recently difficult at best (21,24). The rise of next-generation sequencing and single-molecule real-time (SMRT) technologies tremendously improved the quality of methylome analyses, but it also made it much more accessible (25,26). SMRT sequencing uses a sequencing-by-synthesis approach with fluorescently labeled nucleotides. Pulse width, the signal of nucleotide incorporation, and interpulse duration, the time between two incorporations, allow discrimination between incorporated bases and their methylation status (27).
The purpose of the present study is to understand how DNA methylation regulates gene expression in B. cenocepacia. To this end, a genome-wide methylome analysis was carried out, and genes under DNA methylation regulation were identified.

RESULTS
Identification of B. cenocepacia DNA MTases. All predicted DNA MTase genes in the B. cenocepacia J2315 genome were identified using REBASE (see Table S1 in the supplemental material). DNA MTase genes BCAL3494 and BCAM0992, widely distributed within the genus Burkholderia, were selected for further analysis. Gene BCAL3494, located on the first replicon of B. cenocepacia, is a type III methyltransferase that is part of an RM system, together with a restriction enzyme encoded by the neighboring gene BCAL3493. Gene BCAM0992 is located on the second replicon and apparently does not have any adjacent genes coding for restriction enzymes. Instead, BCAM0992 is part of the trp operon that comprises genes involved in tryptophan biosynthesis and is located between trpB (BCAM0991) and trpA (BCAM0993). It is potentially under transcriptional control of proteins regulating this operon, including the trp repressor as described in E. coli (28). The presence of MTase genes in trp-containing operons has been reported before, as in most enterobacteria the gene of the adenine-specific MTase Dam is located close to trpS (29). Although BCAM0992 was identified as an orphan MTase gene, it strongly resembles MTase genes of type II RM systems, in which the restriction and modification enzymes act separately and are not dependent on each other (30). To investigate the influence of BCAL3494 and BCAM0992 on bacterial physiology, deletion mutants were constructed (Fig. S1). For the other DNA MTase genes in B. cenocepacia J2315 identified with REBASE (Table S1), homologues in different Burkholderia strains could not be found; these genes were not further investigated in the present study.
Phenotype of mutant strains. MTase genes BCAL3494 (including neighboring genes BCAL3488 to -3492) and BCAM0992 were deleted in B. cenocepacia strains J2315 and K56-2 (Fig. S1), and the phenotype of the deletion mutants was investigated in detail. No differences in growth rate during exponential phase were observed between wild-type and mutant strains when cultured in phosphate-buffered mineral medium (Fig. S2). Microscopic analysis showed a different, more clustered biofilm morphology for both BCAL3494 deletion mutants (ΔBCAL3494) compared to wild-type strains, whereas the biofilm structure of the BCAM0992 deletion mutants (ΔBCAM0992) did not differ from wild type (Fig. 1A). Cell aggregation in planktonic cultures was investigated using flow cytometry (Fig. 1B). The degree of aggregation in the BCAL3494 mutant strains was significantly higher (P value J2315, 0.049; P value K56-2, 0.001) than in the corresponding wild-type strains. Also, the ability to form a pellicle, a biofilm-like structure at the air-liquid interface, was investigated (Fig. 1C). Pellicle formation was increased for both ΔBCAL3494 mutants compared to wild-type strains and to ΔBCAM0992 mutants. Complemented mutant strains cΔBCAL3494 and cΔBCAM0992 did not differ significantly from wild type in these experiments, as shown for strain J2315 in Fig. S3.
Motility of all strains was assessed using a swimming motility assay on 0.3% agar plates. After 24 h (strain K56-2) and 32 h (strain J2315), plates were photographed and diameters were measured (Fig. 2). Diameters were significantly smaller for both ΔBCAM0992 mutants compared to wild type (P value J2315, 0.002; P value K56-2, Ͻ0.001). Both ΔBCAL3494 mutants, as well as the complemented mutants, were identical to wild type. We also investigated swarming motility, but no significant differences between the different strains were observed.
Galleria mellonella (wax moth) larvae were used as infection model to assess the virulence of B. cenocepacia J2315 in the absence of methylation. After both 24 h and 48 h, the percentage of survival in the ΔBCAL3494 mutant group was higher than in the other groups. However, the observed difference was not statistically significant (P value, 0.321). After 72 h, all larvae were reported dead (Fig. S4). Although DNA methylation depletion did not lead to a decreased virulence in G. mellonella, the role of adaptive immunity found in higher organisms and the behavior of immune cells toward such MTase mutants with altered biofilm production or swimming motility remains to be investigated.
Effect of the DNA MTase inhibitor sinefungin on methylation-dependent phenotypes. Sinefungin, a structural analog of SAM and known for blocking base methylation in other bacteria such as Streptococcus pneumoniae (31), was used as a DNA MTase inhibitor. The MIC of sinefungin in B. cenocepacia J2315 and K56-2 was determined and was higher than 200 g/ml. Both strains were exposed to sinefungin concentrations below the MIC of sinefungin (50 g/ml) to ensure that any effect observed was not due to growth inhibition by sinefungin, and the effect on biofilm formation, pellicle formation, cell aggregation, and motility was quantified (Fig. 3). Bacteria exposed to sinefungin produced more pellicle mass, showed a higher degree of cell aggregation (P value, 0.003), had a different biofilm morphology, and were less motile (P value, 0.004). These findings suggest that chemically blocking DNA methylation or deleting genes responsible for DNA methylation leads to the same phenotypes in B. cenocepacia J2315 and K56-2.
Methylome analysis. Using SMRT sequencing (PacBio), the complete methylome of B. cenocepacia J2315 and K56-2 was identified. Only data for strain J2315 are reported below, as data for strain K56-2 were highly comparable ( Fig. 4 and 5). Two distinct methylated motifs were identified in the wild-type strain: CACAG and GTWWAC. The CACAG motif was methylated at the fourth position on the forward strand, whereas the GTWWAC motif was methylated at the fifth position on both the forward and reverse strands. Although all CACAG and GTWWAC motifs were methylated in the wild-type strains, methylation of the CACAG motif was absent in the ΔBCAL3494 deletion mutants, and likewise, no methylation of the GTWWAC motif was seen in the ΔBCAM0992 mutants (Table S1). This demonstrates that MTase BCAL3494 recognizes the CACAG motif, while MTase BCAM0992 recognizes the GTWWAC motif.
The location of every methylated CACAG and GTWWAC motif was mapped ( Fig. 4  and 5). In all, 6,834 methylated CACAG motifs and 961 methylated GTWWAC motifs were found, of which the majority was present on the first replicon (CACAG, 45.6%; GTWWAC, 49.9%), followed by the second replicon (CACAG, 42.1%; GTWWAC, 38.9%), the third replicon (CACAG, 10.6%; GTWWAC, 9.0%), and the plasmid (CACAG, 1.7%; GTWWAC, 2.2%). Subsequently, all genes with methylated motifs in their promoter region, here defined as 60 bases upstream of the transcription start site, were identified. Ninety-one promoter regions contained methylated CACAG motifs, and 80 promoter regions contained methylated GTWWAC motifs, with most of the motifs being present on the first replicon ( Fig. 4 and 5). To reveal possible sites of epigenetic regulation, at which DNA methylation is hindered by binding of a regulator, positions of the few unmethylated motifs were analyzed. However, none of these motifs were found in promoter regions. Functional classes of genes found in the data set of genes with methylated promoters include genes involved in intermediary metabolism, regulation, and transport (Table S2).
Virtual Footprint was used to elucidate to which transcription factor (TF) binding sites the discovered methylation motifs CACAG and GTWWAC showed any similarity. Data output of the analysis is listed in Table S3. Sequences that contain methylation motif CACAG were similar to the binding site of E. coli K-12 GlpR, while GTWWACcontaining sequences were similar to binding sites of several other E. coli K-12 TFs, including ArcA, Fis, OxyR, and Fur. Interestingly, the latter two have previously been described as TFs with DNA methylation-blocking ability in Salmonella enterica and E. coli (32,33).
Expression of genes with a methylated promoter. The expression level of genes with methylated promoter regions was determined in wild-type and mutant strains, using qPCR. Expression data for genes with methylated promoter regions are listed in Tables 1 and 2. Volcano plots (Fig. S5) (fold changes plotted against corresponding P values) show that most genes tested were upregulated in the mutants compared to the wild-type strains. Six of these genes were significantly upregulated in mutants of both strain backgrounds: BCAL1515, BCAL2465, and BCAM0820 were upregulated in ΔBCAL3494, whereas genes BCAL0079, BCAL2415, and BCAM1362 were upregulated in ΔBCAM0992. Four additional genes were upregulated in K56-2 mutants only: BCAL0423, BCAM2738, and BCAS0223 were upregulated in ΔBCAL3494, and BCAL1556 was upregulated in ΔBCAM0992. Subsequently, the methylated promoter regions of these genes were analyzed in detail (Fig. 6). In most cases, the methylated motif was in close proximity to the Ϫ10 or Ϫ30/35 element in bacterial promoter regions. As DNA methylation in bacteria can control binding of TFs outside the promoter region (14), we extended the genomic context to look for additional methylation sites upstream of these upregulated genes with a methylated promoter (up to 200 bp), but no additional methylation patterns could be identified.
To confirm that the presence of methylation close to the Ϫ10 or Ϫ30/35 element influences transcription and therefore gene expression in B. cenocepacia, translational The total number of methylated CACAG motifs and methylated CACAG motifs in promoter regions, per replicon (red, replicon 1; green, replicon 2; yellow, replicon 3; gray, plasmid), is shown on the large and small inner circle, respectively. The positions and names of genes with methylated promoter regions are indicated with colored labels (same color code). enhanced green fluorescent protein (eGFP) reporter fusions were constructed and eGFP production was quantified. The eGFP production in strains harboring different plasmids is shown in Fig. 7. As expected, the production of eGFP, driven by the promoters of genes BCAL1515, BCAM0820, and BCAL0079, was significantly (P ϭ 0.001, P ϭ 0.014, and P ϭ 0.002, respectively) increased in the deletion mutant for which an upregulation of these genes was observed using qPCR experiments (Fig. 7).
DNA methylation in the origin of replication. DNA methylation was detected in all origins of replication of B. cenocepacia (Fig. 8). Similar methylation patterns were observed in the origins of the different replicons. A previously discovered 7-mer (CTGTGCA) that can be found in all replication origins (34) contains a CACAG methylation motif on the antisense strand. This motif was also found at the 3= end of almost every DnaA box. These boxes are bound by DnaA proteins, essential for DNA unwinding and chromosome replication initiation (35). Also, the GTWWAC motif was found in proximity to the replication origins; consequently, the origins in B. cenocepacia represent methylation-rich regions. Whereas methylated CACAG motifs were found throughout the origins of replication, the position of the GTWWAC methylation was unique in The total number of methylated GTWWAC motifs and methylated GTWWAC motifs in promoter regions, per replicon (red, replicon 1; green, replicon 2; yellow, replicon 3; gray, plasmid), is shown on the large and small inner circle, respectively. The positions and names of genes with methylated promoter regions are indicated with colored labels (same color code).
all replicons and at least two GTWWAC motifs were found in between two CACAG methylated DnaA boxes. In contrast to the origins of the three larger replicons, the origin of replication of the plasmid contained only one CACAG methylated DnaA box.
To assess the impact of DNA methylation within the replication origins on replication, we visualized the nucleoids in B. cenocepacia wild-type and mutant cells at different time points and calculated the relative nucleoid size (Fig. 9). In both ΔBCAL3494 (49.7%) and ΔBCAM0992 (48.7%), this ratio significantly differed from wild type (33.8%, P value of Ͻ0.001), which suggests that methylation of the replication origins influences DNA compaction. In general, this observation strengthens the proof of replication regulation by DNA methylation. We also investigated differences in cell morphology, but no differences in shape, structure, or size could be observed.

DISCUSSION
Despite the growing knowledge of DNA methylation in prokaryotes (15), the role of DNA MTases in regulating gene expression in B. cenocepacia remains to be revealed. In the present study, we identified two DNA MTases (BCAL3494 and BCAM0992), and mutants in which these genes were deleted showed differences in biofilm formation   (36), the same phenotypic differences were observed. SAM is important in cysteine metabolism (37), and while it cannot be ruled out that addition of sinefungin has an effect beyond DNA methylation, the observation that only differences in biofilm, cell aggregation, and motility were observed (i.e., the same phenotypes affected by MTase deletion) suggests this is not the case. These findings demonstrate that epigenetic control of gene expression by MTases plays an important role in controlling certain phenotypes. Similar results have been reported in Salmonella enterica, where DNA methylation is crucial for optimal pellicle and biofilm production (38). The phenotypic changes in biofilm formation and cell clustering did not correlate with a significant decrease in virulence of B. cenocepacia in the G. mellonella infection model, suggesting DNA methylation is not essential for full virulence in B. cenocepacia.
Methylome analysis showed that mutants in which MTase ΔBCAL3494 or ΔBCAM0992 was inactivated lacked adenine methylation in specific motifs. MTase BCAL3494 was specifically linked to methylation of the CACAG motif, and MTase BCAM0992 was linked to methylation of the GTWWAC motif. This strategy of DNA methylation analysis, in which the methylome of strains lacking MTases is determined, has been used in various bacteria, as it is an effective way to find associations between predicted MTases and genome-wide methylation motifs (39,40). For example, several methylation motifs were identified in Burkholderia pseudomallei, including motifs CA-CAG and GTWWAC (41). Two of the B. pseudomallei MTases (M.BpsI and M.BpsII) are homologous to the B. cenocepacia MTases BCAL3494 and BCAM0992. In Ralstonia  solanacearum, an important plant pathogen that is phylogenetically related to B. cenocepacia, the GTWWAC methylation motif cooccurs with the respective homolog of the BCAM0992 MTase, whereas a BCAL3494 MTase homolog and methylation of CACAG are absent (42). As in B. cenocepacia, the BCAM0992 homolog in R. solanacearum is an orphan DNA MTase. Analysis of cytosine methylation suggests that cytosine is more likely methylated at random instead of at specific motifs and is likely  not having a major regulatory function. Also, GC-rich genomes complicate the search for specific cytosine motifs.
Previous epigenetic research demonstrated that in most cases, there is a negative correlation between methylation in promoters and transcription (43). To uncover the role of DNA methylation in regulation of B. cenocepacia gene expression, all methylated motifs in promoter regions were identified. The data obtained in the present study indicate that gene expression was upregulated in DNA MTase mutants, suggesting that adenine DNA methylation in B. cenocepacia affects gene expression by a mechanism inhibiting transcription. In both prokaryotes and eukaryotes, adenine and cytosine methylation are involved in blocking (or enhancing) the binding of RNA polymerase to DNA (15,21,44), and especially methylation near the Ϫ10 and Ϫ30/35 elements in the promoter region may affect RNA polymerase binding (45). We found that also in B. cenocepacia, methylated motifs (CACAG and GTWWAC) are found close to or in these elements.
BCAM0820, upregulated in the J2315 and K56-2 ΔBCAL3494 mutant, is a twocomponent response regulator, the first gene of an operon homologous to the Wsp chemosensory system involved in biofilm formation in Pseudomonas aeruginosa (46). BCAM0820 is homologous to WspR but lacks the diguanylate cyclase domain. During an experimental evolution study in which B. cenocepacia HI2424 biofilms were grown on beads, mutations within the wsp gene cluster occurred in different clones; these were associated with increased pellicle formation and increased biofilm formation on beads. This demonstrates that the Wsp cluster is involved in pellicle formation in B. cenocepacia (47,48), and the upregulation of BCAM0820 could explain the differences in pellicle and biofilm formation between the wild-type strains and the ΔBCAL3494 deletion mutants observed in the present study. Interestingly, BCAL1515, encoding 2-oxoglutarate dehydrogenase (SucA) and upregulated in ΔBCAL3494, also acquired mutations in the course of the experimental evolution study (48), but the role of this gene in biofilm formation has not been further explored. BCAL0079, upregulated in the ΔBCAM0992 mutants, is annotated as a DNA helicase gene (rep). Besides unwinding DNA during DNA replication, Rep plays a role in swimming motility in E. coli (49). The reduced motility observed in the ΔBCAM0992 mutants suggests that Rep may also affect motility in B. cenocepacia, although this remains to be confirmed.
Measurement of eGFP production in translational fusion mutants revealed that mutants with constructs containing the BCAL1515, BCAM0820, or BCAL0079 promoter showed a significant increase in eGFP production compared to wild type, thereby supporting our hypothesis of gene expression regulation by DNA methylation. In silico analyses predict that sequences containing methylation motifs are similar to binding sites of TF in E. coli K-12, and it is plausible that these sequences are also part of TF binding sites in B. cenocepacia, allowing us to propose a possible mechanism of gene expression regulation (Fig. 10). TFs that bind close to the Ϫ10 and Ϫ35 region often act as transcriptional repressors (50). Therefore, a methylated promoter region could promote binding of a repressor (51) and sterically hinder RNA polymerase (OFF state), whereas an absence of methylation could lead to binding of RNA polymerase and initiation of transcription (ON state).
The role of DNA methylation in prokaryotes in multifaceted. Besides gene expression regulation and a role in DNA mismatch repair in Gram-positive bacteria (52), DNA methylation has also been implicated in the coordination of replication initiation. Results of the present study seem to confirm this, as the rep gene, necessary for replication, is under DNA-methylation-mediated epigenetic control. In E. coli, GATC motifs, omnipresent in the replication origin, are prone to adenine methylation. The motifs are found within DnaA boxes, essential for binding of the DnaA protein and initiation of replication. The methylation state of each of these GATC motifs changes the affinity of DnaA and sequestering-protein SeqA for the DnaA box. Immediately after replication, GATC motifs are hemimethylated, which leads to sequestration of the DnaA boxes by SeqA and prevents the reinitiation of DNA replication (21). The occurrence of methylated motifs in the vicinity of the origins of replication of the four replicons in B.
cenocepacia was studied to check for a link between DNA methylation and coordination of the replication process. An enrichment of the CACAG motif was observed in the origin of replication of all replicons. The motif was part of a bigger sequence that has previously been reported as a recurring 7-mer (34), without known function. In addition, the origin of replication of the different replicons showed high similarities in methylation patterns, and both mutants appeared to have a larger nucleoid size with respect to the total cell size compared to wild type, raising the possibility of replication coordination by DNA methylation. The importance of adenine methylation in DNA replication has been studied in other bacteria besides E. coli, and for example in Vibrio cholerae adenine methylation is also frequently found in the origin of replication of both chromosomes (53), just like in B. cenocepacia. Although the lack of differences in growth between B. cenocepacia wild type and MTase mutants during exponential phase suggests the impact of DNA methylation on replication is limited, the K56-2 ΔBCAL3494 mutant did reach maximal optical density in stationary phase a bit earlier than the wild type (see Fig. S2 in the supplemental material). This is in line with the observation that the promoter of dnaA (BCAL0423) was methylated at multiple positions and under transcriptional control of the BCAL3494 MTase (Fig. 6). Combined, these data suggest that DNA methylation contributes to coordinating the replication process in B. cenocepacia K56-2. In the slower-growing J2315 strain, such difference in maximal optical density was not observed, and in this strain, expression of dnaA was found not to be transcriptionally controlled by BCAL3494, suggesting that coordination of replication is strain dependent.
In conclusion, we have demonstrated that DNA methylation plays a role in regulation of gene expression in B. cenocepacia. DNA MTases BCAL3494 and BCAM0992 are essential for methylation of the B. cenocepacia genome and are responsible for methylation of base motifs CACAG and GTWWAC, respectively. The absence of methylation resulted in an upregulation of certain genes with methylatable promoter regions, including BCAM0820 and BCAL0079, genes of which their function can be linked to the observed biofilm-and motility-affected phenotypes. Finally, recurrent methylation patterns were detected in all origins of replication, which suggests an additional role of DNA methylation in replication regulation.

MATERIALS AND METHODS
Strains and culture conditions. All strains and plasmids used in this study are listed in Table 3  ). LB medium (Luria-Bertani medium with 5 g/liter NaCl; Sigma-Aldrich) was used for maintenance of E. coli strains and during specific stages of the gene deletion procedure (see below) where antibiotic selection with tetracycline (250 g/ml) (Sigma-Aldrich) was desired. Prior to phenotypic experiments, liquid overnight cultures were grown in a shaker incubator (100 rpm) at 37°C. Selection of DNA MTase genes-in silico. The REBASE Genome database was used to allocate all known DNA MTase genes in the B. cenocepacia J2315 and K56-2 genomes (54). The Artemis Genome Browser and Annotation Tool (Sanger) allowed visualization of the genomic context of these genes (55). NCBI BLAST was used to screen for conservation of the genes within the Burkholderia genus using default search parameters (56) (search mode, BLASTn; E cutoff value, Ͻ1EϪ5).
Construction of deletion mutants. All primers used for construction and complementation of the deletion mutants are listed in Table S4 in the supplemental material. The procedure is an adapted allelic replacement approach, using a suicide plasmid with an SceI endonuclease recognition site (57,58). The suicide plasmid, containing DNA fragments of regions flanking the target gene, is integrated into the B. cenocepacia genome by homologous recombination. Introducing a second plasmid that carries SceI endonuclease genes into B. cenocepacia results in a lethal genomic strand break. Another homologous recombination event allows the bacteria to repair the break with a 50% chance of resulting in a gene deletion. Deletion mutants ΔBCAL3494 and ΔBCAM0992 were constructed in both B. cenocepacia J2315 and K56-2.
BCAL3494 was deleted together with neighboring gene BCAL3493, as well as BCAL3488 to BCAL3492 (encoding hypothetical proteins). Targeting BCAL3494 alone was not feasible because regions flanking BCAL3494 contain multiple recognition sites for endonucleases used during construction of the deletion mutants, and digestion of these regions would be inevitable (Fig. S1).
E. coli One Shot PIR2 cells (Thermo Fisher), expressing pir, were used for transformation, replication, and maintenance of the suicide plasmid during construction of deletion mutants. Thawed cells were immediately exposed to a heat shock transformation procedure, after which they were transferred to SOC (Super Optimal broth with Catabolite repression) medium for recovery. For plasmid selection, the phosphate-buffered mineral medium was supplemented with one or more of the following antibiotics: trimethoprim (Ludeco; 50 g/ml for initial screening in E. coli, 200 g/ml when plasmid is introduced in B. cenocepacia), chloramphenicol (400 g/ml), gentamicin (50 g/ml), kanamycin (50 g/ml), and ampicillin (200 g/ml) (Sigma-Aldrich).
Construction of plasmids for complementation. To ensure that phenotypes were solely caused by the deletion of DNA MTases, deletion mutants were complemented. The primers used for construction of plasmids used for complementation are listed in Table S4. Plasmids pJH2 and pSCrhaB2 were used for complementation of ΔBCAL3494 (cΔBCAL3494) and ΔBCAM0992 (cΔBCAM0992), respectively. The genomic sequences of the DNA MTase genes were PCR amplified and subsequently cloned into the plasmids. BCAL3494 was amplified, including its own regulatory region (approximately 250 nucleotides upstream of the transcription start site), into pJH2, which does not have a promoter associated with its multiple cloning site (59). BCAM0992 does not have its own upstream promoter; therefore, it was cloned into pSCrhaB2, which contains a rhamnose-inducible promoter (60). Complemented mutant strains were subjected to the same phenotypic tests as the deletion mutants and wild-type B. cenocepacia. For strain cΔBCAM0992, the phosphate-buffered mineral medium was supplemented with 0.2% rhamnose.
Biofilm and clustering experiments. Biofilms were grown in plastic U-shaped 96-well microtiter plates in phosphate-buffered medium at 37°C, starting from 200-l/well planktonic overnight cultures with an optical density (OD) of 0.05 (590 nm). After 4 h static incubation, all wells were rinsed with physiological saline (PS; 0.9% NaCl in water), thereby removing all unattached planktonic cells. Wells were refilled with 200 l medium and incubated for an additional 20 h. Where appropriate, biofilms were stained with LIVE/DEAD (SYTO9/propidium iodide; Invitrogen) to visualize the bacteria and distinguish live and dead cells (61). Pellicle formation was determined in glass tubes. Cultures were grown statically for 24 h, after which adhering pellicles were stained and quantified with crystal violet (62). Cell clustering, already shown to be correlated with pellicle formation, was determined with flow cytometry (Attune NxT flow cytometer; Thermo Fisher) (63). Forward scatter (FSC), a value for particle size, and side scatter (SSC), a value for particle complexity, were measured for each particle present in the bacterial suspension and visualized in scatterplots. After analysis of these graphs, the main cell population was gated (gate ranging from approximately 10 3 to 10 5 for both FSC and SSC), and detected events larger and more complex than the gate were considered clustered (Fig. S6).
Motility experiments. Petri dishes containing phosphate-buffered mineral medium with agar concentrations of 0.3% and 0.5% were used for assessment of swimming and swarming motility, respectively. One microliter of cultures with an OD of 0.1 was spotted on the agar plates. Diameters were measured after 24 h (strain K56-2) or 32 h (strain J2315).
DNA MTase inhibition with sinefungin. A stock solution of the DNA MTase inhibitor sinefungin (Sigma-Aldrich) was prepared (10 mg/ml) (31), aliquoted, and immediately frozen at Ϫ20°C to prevent degradation. Cells were grown for 24 h in sinefungin-supplemented medium (50 g/ml) and used as inoculum for an overnight culture, also in sinefungin-supplemented medium. This allowed the DNA MTase inhibitor to have an effect during several growth cycles. Then, biofilm formation and motility of sinefungin-treated cells were assessed as described above in medium supplemented with 50 g/ml sinefungin.
Virulence assay. Galleria mellonella larvae (TruLarv; BioSystems Technology) were injected in the left proleg with 10 l bacterial suspension (set to an OD of 0.2 and diluted 100-fold before injection). A physiological saline solution was used as control. After infection, the larvae were incubated at 37°C, upon which viability was checked every 24 h. Eight larvae were included in each group.
Genomic DNA extraction. Prior to DNA extraction, planktonic strains were grown overnight in a shaker incubator (100 rpm) at 37°C. Biofilm cells were grown as described above. Next, cells were harvested and genomic DNA (gDNA) was extracted using the Wizard Genomic DNA purification kit (Promega). Quantification was performed with a BioDrop LITE (BioDrop) spectrophotometer.
SMRT sequencing. To determine the methylome of B. cenocepacia, gDNA extracts were analyzed with single-molecule real-time (SMRT) sequencing technology. gDNA samples of both wild-type and mutant strains were run on a Pacific Biosciences Sequel System (250ϫ coverage) according to the manufacturer's guidelines. Library preparations were multiplexed as data output of approximately 2 Gb per genome was expected, and a single SMRT Sequel cell provides up to 6 Gb data. Initial data output was processed with SMRT Link software (Pacific Biosciences). Identification of the modified bases and analysis of the methylated motifs were performed with the Base Modification and Motif Analysis application (SMRT Link v6.0; Pacific Biosciences). In-depth data analysis was performed with CLC Workbench Genomics (v11.0.1; Qiagen). Differential analysis between wild type and mutants was performed to identify methylation motifs specifically associated with certain DNA MTases. Previously predicted promoter regions and transcription start sites of B. cenocepacia were used to determine the methylation profile of regulatory regions (64). Virtual Footprint software (promoter analysis mode, default search parameters) was used to assess similarity of the methylation motifs to known TF binding sites (65).
qPCR. To evaluate the impact of DNA methylation in promoter regions on gene expression, qPCR was performed on all genes that had a methylated promoter region in wild-type B. cenocepacia but an absence of methylation in the promoter region in one of the deletion mutants. All hypothetical genes and genes with unknown function, as well as genes with low innate expression level, were excluded from testing. The primers used for qPCR are listed in Table S4. First, all strains were grown to an OD of 0.6 in phosphate-buffered medium, after which they were pelleted by centrifugation and frozen at Ϫ80°C. Next, RNA was extracted using the RiboPure bacterial extraction kit (Invitrogen), followed by a DNase treatment to remove trace quantities of gDNA. Quantification and measurement of RNA purity of the extracts were performed with a BioDrop LITE (BioDrop). Subsequently, cDNA was synthesized, using 500 ng RNA per reaction, with a reverse transcriptase kit (high-capacity cDNA RT kit; Applied Biosystems). Per qPCR mixture, 2 l template cDNA was mixed with 10 l GoTaq qPCR master mix, 0.6 l qPCR primer mix (10 g/ml), and 7.4 l nuclease-free water according to the GoTaq qPCR master mix (Promega) protocol. Samples were run on a CFX96 Real-Time System C1000 Thermal Cycler (Bio-Rad), and output data were processed with Bio-Rad CFX Manager 3.1 software. The baseline threshold was set to a defined 100 relative fluorescence units. Obtained quantification cycle (C q ) values were normalized to reference gene rpoD (BCAM0918) of which the expression was stable across all samples; differences from wild type were calculated (ΔΔC q ) and log-transformed. Volcano plots were used to plot the negative logarithm of statistical P values against log 2 fold changes (Fig. S5).
Construction of translational eGFP reporter fusions and measurement of eGFP production. Genes with methylated promoter regions that showed a significant upregulation of gene expression in one of both mutant strains were selected for eGFP experiments. Translational eGFP reporter fusion plasmids were constructed by cloning the regulatory regions of the genes, comprising 60 to 390 nucleotides upstream of the transcription start site, into vector pJH2. The insert is cloned right in front of the eGFP gene and contains an ATG start codon at the 3= end, in frame with the codon sequence of the gene. All primers used for amplification of the regulatory regions and screening of pJH2 with correct insert length are listed in Table S4. The plasmids were transferred to B. cenocepacia J2315 and K56-2 by triparental mating. Exconjugants were grown on selective plates (LB medium supplemented with 200 g/ml chloramphenicol and 50 g/ml gentamicin) and PCR screened to confirm the presence of the insert. Constructs carrying genes BCAL2415, BCAL2465, and BCAM1362 repeatedly failed to be transferred to B. cenocepacia and were not included in further experiments. Fluorescent signals of eGFP production in wild-type and mutant strains were measured by flow cytometry (Attune NxT flow cytometer; Thermo Fisher) (59).
Nucleoid staining and surface quantification. To stain the nucleoids of B. cenocepacia, 4=,6diamidino-2-phenylindole (DAPI) (Thermo Fisher) was used. In brief, cells were grown for various times, heat fixed on glass slides, stained with a DAPI solution (3 g/ml) for 5 min, and visualized with fluorescence microscopy (EVOS FL Auto Imaging System; Thermo Fisher). To calculate nucleoid and total cell surfaces, we used the thresholding function of the Java-based processing application ImageJ (LOCI, University of Wisconsin).
Data analysis and statistics. Statistical analysis was performed using SPSS Statistics v. 25 software. All tests and experiments were run in triplicate unless otherwise mentioned. Normality of data was verified with a Shapiro-Wilk test. To check for significant differences between data, normally distributed data were subjected to a t test or one-way analysis of variance (ANOVA) test, not normally distributed data to a nonparametric Mann-Whitney U-test. A log rank (Mantel-Cox) test was used to assess G. mellonella data. Resulting P values smaller than 0.05 were reported as statistically significant.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.