N4-Cytosine DNA Methylation Is Involved in the Maintenance of Genomic Stability in Deinococcus radiodurans

DNA methylation serves as a vital component of restriction-modification (R-M) systems in bacteria, where it plays a crucial role in defense against foreign DNA. Recent studies revealed that DNA methylation has a global impact on gene expression. Deinococcus radiodurans, an ideal model organism for studying DNA repair and genomic stability, possesses unparalleled resistance to DNA-damaging agents such as irradiation and strong oxidation. However, details on the methylome of this bacterium remain unclear. Here, we demonstrate that N4-cytosine is the major methylated form (4mC) in D. radiodurans. A novel methylated motif, “C4mCGCGG” was identified that was fully attributed to M.DraR1 methyltransferase. M.DraR1 can specifically bind and methylate the second cytosine at N4 atom of “CCGCGG” motif, preventing its digestion by a cognate restriction endonuclease. Cells deficient in 4mC modification displayed higher spontaneous rifampin mutation frequency and enhanced DNA recombination and transformation efficiency. And genes involved in the maintenance of genomic stability were differentially expressed in conjunction with the loss of M.DraR1. This study provides evidence that N4-cytosine DNA methylation contributes to genomic stability of D. radiodurans and lays the foundation for further research on the mechanisms of epigenetic regulation by R-M systems in bacteria.


INTRODUCTION
In addition to the canonical four-base hereditary information, numerous chemical modifications occur in DNA nucleotides, vastly expanding its structural complexity and information depth (Korlach and Turner, 2012). Among them, DNA methylation is a ubiquitous and pivotal epigenetic phenomenon that occurs in the genomes of all organisms. In bacteria, DNA methylation serves as a vital component of restriction-modification (R-M) systems, which plays a crucial role in defense against foreign DNA (Loenen et al., 2014;Pingoud et al., 2014;Rao et al., 2014;Pleska et al., 2016). It can also modulate the interactions of DNA-binding proteins with their cognate sites and is involved in chromosome replication, correction of DNA mismatches, cell cycle-coupled transcription, and formation of epigenetic lineages (Adhikari and Curtis, 2016;Casadesus, 2016). Recent studies demonstrated that DNA methylation had a global impact on gene expression (Casselli et al., 2018;Kumar et al., 2018).
N 6 -methyladenine (6mA), N 4 -methylcytosine (4mC), and C 5methylcytosine (5mC) are the major methylated sites in bacterial genomes, which are generated by methyltransferases (MTases) transferring a methyl group from S-adenosyl-l-methionine (SAM) to certain target bases. To date, the overwhelming majority of sequenced bacterial genomes have been found to encode MTases and are subject to DNA methylation Blow et al., 2016). However, the frequency and distribution of methylated bases, as well as the precise sequence targets and biological roles of most MTases remain unknown, largely due to a lack of efficient analytical methods for the detection of these modified bases. Mass spectrometry has long been established as a convenient, rapid and effective technique for qualitative and quantitative analysis of methylated bases (Greer et al., 2015;Schiffers et al., 2017). However, it is unable to achieve genome-wide mapping of target sequences. Fortunately, new cutting-edge methods provide great promise in addressing this issue. Single-molecule, real-time sequencing (SMRT-seq) and nanopore-based sequencing have expanded the epigenetics research toolkit, offering novel methods for the genome-wide, systematic detection of the three forms of methylated DNA and unprecedented opportunities to explore the spectacular landscape of bacterial epigenomes Fang et al., 2012;Rand et al., 2017).
Deinococcus radiodurans is a robust bacterium with an unparalleled DNA repair system. Experimental evidence has shown that D. radiodurans has a similar rate of doublestrand break (DSB) formation as Escherichia coli under identical conditions (Gerard et al., 2001). However, after receiving a dose of 5,000 gray of γ-irradiation, which introduces approximately 1,600 DSBs per cell, D. radiodurans maintains its viability and is able to rebuild its shattered genome within hours (Cox and Battista, 2005). Several mechanisms explain the resistance of D. radiodurans to extreme radiation and oxidative stress, including strong self-regulatory systems, cell defense systems and efficient DNA repair mechanisms. In contrast to the longheld assumption that D. radiodurans DNA was methyl-deficient due to the unidentified methylated bases or undetectable DNA MTase activity (Ferrandi et al., 2019), the genomic DNA (gDNA) of D. radiodurans was reported to contain 6mA and 5mC modifications (Prasad et al., 2005;Patil et al., 2017). The hypomethylation of 6mA leads to dysregulation of metabolic pathways (Shaiwale et al., 2015), whereas 5mC could regulate protein homeostasis in this bacterium (Patil et al., 2017). However, recent experiments using dot-blot and enzyme-linked apta-sorbent assays revealed the absence of 5mC modifications in the gDNA of D. radiodurans (Ferrandi et al., 2019).
In the present study, we utilized SMRT-seq combined with sensitive ultra-high-performance liquid chromatography coupled with triple-quadrupole tandem mass spectrometry (UHPLC-QQQ-MS/MS) to analyze DNA methylation patterns in D. radiodurans. A novel 4mC modification motif was identified and fully attributable to the presence of the M.DraR1 gene.
M.DraR1 could methylate the target motif and prevent it from being digested by a cognate restriction endonuclease (REase) both in vivo and in vitro. Cells lacking M.DraR1 were assessed on their varied phenotypes and gene expression patterns. A large number of differentially expressed genes (DEGs) were observed in conjunction with the loss of MTase, including genes involved in DNA repair, transcription and nucleotide processing, suggesting that 4mC methylation is involved in the regulation of gene expression and genomic stability in D. radiodurans.

RESULTS
N 4 -Cytosine Is the Major Methylated Form in D. radiodurans Deinococcus radiodurans R1 strain contains five type II and two type IV R-M systems (Supplementary Figure S1 and Supplementary Table S3). A fused R-M enzyme encoded by the open reading frame (ORF) ORF2230P, DraRI was previously characterized as a member of the MmeI family recognizing the sequence CAAGN 6m AC and modifying the penultimate adenine (Morgan et al., 2009). ORF14075P encoding MTase belongs to the YeeA (COG1002) methylase, which shares 32% similarity with the protein DrdIV that recognizes and methylates the sequence TACG 6m AC (Morgan et al., 2009). The MTase encoded by ORF15360P contained the typical motifs associated with a 6mA MTase, but its recognition sequence could not be predicted because it lacked the similarity to other MTases with identified specific substrate. ORF16000P encodes a putative MTase, designated as M.DraR1, which has been reported to be an atypical SAM-dependent 5mC MTase for recognizing and methylating random 5m CpN sites (Patil et al., 2017).
Genome-wide analysis of all template positions revealed a cluster of cytosine residues that were clearly separated from the background (Figures 1A,B). Motif analysis identified two MTase motifs, "C 4m CGCGG" and "DNCGGTGTGNB" with an unknown modification at the fifth base ( Figure 1C and Table 1). And no 6mA modification motifs were identified in the gDNA of D. radiodurans R1. In addition, the abundance of methylated bases was quantified by UHPLC-QQQ-MS/MS, and results showed that the gDNA of D. radiodurans contains approximately 1.3 molar ratio of 4mC/C and 0.4 of 6mA/A (Figure 2 and Supplementary Figure S2), which corresponds to the results of SMRT-seq analysis. No 5mC modification was detected in this bacterium by either method. Taken together, our study suggests that N 4 -cytosine is the major methylated residue in D. radiodurans.
M.DraR1 Is a Novel α-Class N 4 -Cytosine Methyltransferase M.DraR1, located at the small plasmid of genome (A2G07_16000 or DR_C0020) in D. radiodurans, is predicted to be a typical α-class DNA MTases containing a canonical N-terminal SAMbinding motif ("FxGxG") followed by a centrally located target recognition domain and a C-terminal catalytic motif ("SPPY"). Residues in M.DraR1, including V44, D63, P64, T71, G79, K98,  T255, N261, G394, and R404, are also moderately conserved in its homologs (Supplementary Figure S3), indicating that these residues might play important roles in MTase activities. Sequence analysis revealed that M.DraR1 is an N 6 -adenine/N 4 -cytosine methylase. However, the target sequence context, the type and exact site of methylation are still unclear.
To rigorously assess the ability of M.DraR1 transferring a methyl group to either N 4 -cytosine or N 6 -adenine, we FIGURE 2 | Quantification of 6mA, 4mC, and 5mC in the gDNA of D. radiodurans strains using UHPLC-QQQ-MS/MS. The mole ratios are shown. Each column represents the mean and SD of three biological replicates per group. Data with capital letters are significantly different at p < 0.01. Compared with DraR1 WT, * p < 0.05, * * p < 0.01, * * * p < 0.001. generated an M.DraR1 gene knockout strain ( Supplementary  Figures S4A,B). The other three DNA MTase genes were also deleted ( Supplementary Figures S5A-C). Compared to wild type D. radiodurans, deletion of the M.DraR1 gene drastically decreased the molar ratios of 4mC/C, but had no effect on the abundance of 6mA (Figure 2), suggesting that M.DraR1 is an N 4 -cytosine DNA MTase but not the N 6 -adenine MTase. Kinetic signatures obtained from SMRT-seq showed the complete loss of the second 4mC methylation site in the sequence "CCGCGG" of M.DraR1 mutant ( Table 1). The other three DNA MTases were confirmed to be involved in 6mA modification. These results suggest that M.DraR1 specifically recognizes the target sequence "CCGCGG" and methylates the N 4 atom of the second cytosine in D. radiodurans, indicating it is a novel α-class N 4cytosine MTase.
M.DraR1 Specifically Methylates and Protects the "CCGCGG" Motif From Cleavage by a Cognate REase M.DraR1 protein was purified by stepwise affinity column purification as described in "Materials and Methods" section ( Supplementary Figures S6A-C). And the purified recombinant protein was confirmed to be the product of ORF16000 (Supplementary Figures S6D,E and Supplementary Table S4). M.DraR1 exists mainly in its monomeric state in solution during the size exclusion column (Supplementary Figures S6F,G), which is consistent with that of previously reported DNA MTases (Ma et al., 2016). To test the DNA binding capacity of the M.DraR1 protein, different oligonucleotides were generated by annealing corresponding single strands (Supplementary Table S2). As shown in Figure 3, M.DraR1 was able to bind oligonucleotides containing "CCGCGG" in the presence of SAM  Table S2, from Supplementary Tables S2-S7 paired oligomers). Another blunt-ended dsDNA for "CCGCGG" was analyzed as well (last lane, S8 oligonucleotide cassette). The 100 µM double strand cassettes were used for each reaction. All the reactions were performed in three independent biological replicates.
( Figure 3A), and the binding was saturated by increasing the molar ratio of M.DraR1 to oligonucleotides ( Figure 3A). Furthermore, M.DraR1 was unable to bind oligonucleotides containing any mutations (Figure 3B), which indicates the high specificity of M.DraR1.
To further explore the specificity of M.DraR1 and assess its ability to effectively inhibit REases, restriction analyses were carried out using the pRRS-M.DraR1 vector ( Figure 4A). Methylated pRRS-M.DraR1 plasmid could be digested by HindIII ("AA ↓ GCTT") but not by SacII (a relevant cognate REase), both of them are sensitive to cytosine modifications ( Figure 4B). These results indicate that M.DraR1 cannot methylate "AAGCTT" sites at C residues in vivo but does specifically methylate "CCGCGG" motifs. Other four restriction enzymes with defined specific sequences (CAG ↓ CTG for PvuII, GCG ↓ C for HhaI, C ↓ CGG for HpaII, and GG ↓ CC for HaeIII), which are resistant to cleavage by cytosine methylation at CpN sites, were also analyzed. pRRS-M.DraR1 and its PCR products exhibited similar gel patterns to that of pRRS null vector ( Figure 4B), which suggests that M.DraR1 cannot randomly methylate CpN sites in vivo. MTase activity of M.DraR1 in vitro showed that those DNA substrates were resistant to cleavage by SacII but completely digested by other enzymes (Figure 5 and Supplementary Figure S7). These results also suggested that recombinant M.DraR1 could specifically methylate "CCGCGG" motifs in vitro and protect them from cleavage.

Cells Devoid of M.DraR1 Causes Genomic Instability
Deletion of the M.DraR1 gene did not significantly affect cell growth under optimal conditions ( Figure 6A). The D. radiodurans M.DraR1 mutant showed increased spontaneous rifampin mutation frequency ( Figure 6B). Compared to the wild type, the mutant showed a ∼ninefold increase significantly (p = 0.002), suggesting that M.DraR1 mediated DNA 4mC methylation might be involved in the genomic stability in D. radiodurans. Furthermore, cells lacking M.DraR1 led to a >100-fold increase in transformation frequency ( Table 2). And the relative intermolecular recombination efficiency of the mutant strain was approximately ∼39-fold higher than that of the wild type strain ( Table 2). Compared with the pRADK vector, the modified plasmid pRADKm caused reduced transformability, while the use of methylated M.pRADKm plasmid with M.DraR1 lead to an efficient improvement of DNA transformation frequencies in D. radiodurrans strains ( Table 3). These results suggest that the spontaneous rifampin mutation frequency, recombination capacity and transformability of D. radiodurans might be regulated by 4mC modification. Thus, we proposed that 4mC methylation could act as a methylation barrier against foreign DNA, contributing to genomic stability.

Absence of N 4 -Cytosine DNA Methylation Leads to Expression Changes of Proteins Involved in DNA Damage Response
RNA-seq was used for differential expression gene analysis between wild-type and M.DraR1 strains. A total of 1,158 genes (766 upregulated genes and 392 downregulated genes) were differentially expressed in M.DraR1 strain ( Figure 7A). Hierarchical clustering was performed to determine the expression patterns for all DEGs and samples were clustered into two groups ( Figure 7B). The correlation coefficient of FPKM was >90.1% between two independent biological replicates of each group (Figure 7C), indicating that the replicates were highly correlated and credible. These results were also validated by RT-qPCR ( Figure 7D).
Strikingly, gene ontology (GO) enrichment analysis revealed that upregulated unigenes were significantly enriched in twelve subcategories ( Figure 8A). The majority of unigenes in the biological process category were associated with macromolecule metabolic processes. Interestingly, a number of unigenes were also associated with DNA recombination. Within the molecular function category, most unigenes were assigned to nucleic acid, DNA binding and hydrolase, and nucleoside triphosphatase activity. The KEGG enrichment analytical scatter diagram showed that the major pathways affected in the deletion strain belonged to ABC transport, two-component systems, DNA replication, homologous recombination, and protein export pathways ( Figure 8B). Two-component systems and homologous recombination are the most reported pathways involved in the robust DNA repair capacity of D. radiodurans. In contrast, GO and KEGG pathway enrichment analysis showed that downregulated unigenes were mainly involved in energy or small-molecule metabolic processes, such as the tricarboxylic acid cycle, glycolysis, microbial metabolism in diverse environments, and sulfur metabolism (Supplementary Figure S8).
Interestingly, the majority of DNA damage response genes, such as ddrA, ddrB, ddrO, ddrG, ddrJ, ddrK, ddrD, and pprA, were significantly upregulated ( Table 4). Several genes involved in DNA recombinational repair pathway including recA, recO, rexD, and ruvB were also significantly upregulated. Putative competence genes, including those encoding prepilin peptidase, CRP and CinA, were upregulated in the mutant strain. This gene  expression pattern of the M.DraR1 strain resembled that of D. radiodurans during recovery after exposure to acute radiation (Liu et al., 2003), indicating that the absence of N 4 -cytosine DNA methylation may result in intracellular stress.

DISCUSSION
The complete genome of D. radiodurans contains two circular chromosomes and two circular plasmids with a GC content >66.3% (Hua and Hua, 2016). Despite that the gDNA of D. radiodurans was previously assumed to be methyl-deficient (Ferrandi et al., 2019), Patil et al. detected 6mA and 5mC bases in its gDNA (Prasad et al., 2005;Patil et al., 2017). However, our SMRT-seq and UHPLC-QQQ-MS/MS results suggest that 4mC is the dominant methylated base in D. radiodurans, while the abundance of 6mA is significantly lower and 5mC is undetectable. These results are consistent with those of Ferrandi et al. (2019), who failed to identify any 5mC methylated bases in gDNA of this bacterium. Thus, the N 4 -cytosine DNA methylation appears to be the epigenetic marker in D. radiodurans.
The comparative analysis suggested that MTases encoded by ORF2230P, ORF14075P and ORF15360P are involved in 6mA modification in D. radiodurans, whereas only the MTase encoded by ORF16000P was implicated in 4mC methylation. And the protein encoded by ORF10260P (DR_0643) was annotated as an RsmD-like RNA methyltransferase, which was once proposed to be a 6mA DNA MTase (Shaiwale et al., 2015). As a member of the MmeI family, DraRI in Deinococcus was reported to possess no DNA methyltransferase activity in vivo, but displayed endonuclease activity once disrupted ORF was corrected (Morgan et al., 2009). Levels of 6mA in our study were undetectable by UHPLC-QQQ-MS/MS after deletion of this gene, indicating that DraRI has MTase activity.
The S/DPPY/F motifs are the hallmark of N 4 -cytosine MTase active sites (Malone et al., 1995;Bujnicki and Radlinska, 1999). In contrast, an invariant Pro-Cys dipeptide (PC motif) is present in the active site of 5mC DNA MTases (Kumar et al., 1994). M.DraR1 contains SPPY and FxGxG motifs and but  no other conserved motifs such as PC, ENV, QRR, RE or GN, suggesting that M.DraR1 is not likely a 5mC MTase.
In the current study, deletion of M.DraR1 drastically reduced 4mC content and resulted in complete loss of the modified "C 4m CGCGG" motif, suggesting that M.DraR1 is an N 4 -cytosine MTase. Electromobility shift assay (EMSA) also confirmed that M.DraR1 can rigorously bind to DNA substrates containing the sequence "CCGCGG" in the presence of SAM. Restriction digestion pattern analysis in vivo and in vitro provided further evidence that M.DraR1 is able to protect DNA substrates from SacII digestion. Taken together, these results indicated that M.DraR1 is an α-class N 4 -cytosine MTase, which methylates the N 4 atom rather than C 5 position of the second cytosine in "CCGCGG" motifs (Roberts et al., 2015). Despite the low frequency compared to that of 6mA, 4mC is a prevalent modified base in thermophilic bacterial DNA and also occurs in many mesophilic bacteria (Ehrlich et al., 1985). A recent study revealed that N 4 -cytosine DNA methylation in Helicobacter pylori was involved in transcriptional regulation and pathogenesis, suggesting that 4mC is not only a component of R-M systems but also serves as a global epigenetic regulator (Morgan et al., 2016;Kumar et al., 2018). In our study, 4mC plays a crucial role in the maintenance of genomic stability in D. radiodurans, as higher spontaneous rifampin mutation frequency, enhanced recombination and transformation efficiencies were observed in N 4 -cytosine methylation deficient strains.
Our observations raise the question of why 4mC modifications are the epigenetic marker in D. radiodurans. 4mC has been observed more often in thermophilic bacteria than in nonthermophilic bacteria (Ehrlich et al., 1985;Beaulaurier et al., 2019). This may be related to the fact that N 4 -methylated cytosine is more resistant to heat-induced deamination, which can introduce disadvantageous T:G mismatches into the genome compared with C 5 -methylated cytosine (Ehrlich et al., 1985). Therefore, we speculate that D. radiodurans inherited 4mC methylation from the shared ancestor of Deinococcus-Thermus. We also noted that some other Deinococcus species also possess 4mC and have few 5mC residues in their gDNA (data are not shown), which is consistent with the findings of Blow et al. (2016) who also identified 4mC but no 5mC modifications in three other Deinococcus-Thermus species. Another interesting possibility is that this bacterium may have acquired the MTase gene via horizontal gene transfer since M.DraR1 shows the highest sequence homology with M.Rca13941ORF67P (43% identity and 58% positivity), a type II N 4 -cytosine DNA MTase, identified from Roseiflexus castenholzii DSM13941, were found to recognize the "CCGCGG" motif as well (Roberts et al., 2015).
Natural transformation, the stable uptake and heritable integration of extracellular DNA, plays a major role in genomic diversity, ecological diversification and functional regulation (Thomas and Nielsen, 2005;Johnsborg et al., 2007). In D. radiodurans, the absence of m4C affected the natural transformation pathway through epigenetic regulation of gene expression. This is consistent with our findings that several homologous competence genes, including the genes encoding prepilin peptidase, CRP and CinA were upregulated and the gene encoding extracellular nuclease downregulated in the M.DraR1 strain. Furthermore, uptake of pRADKm DNA was reduced, whereas M.pRADKm was enhanced in both D. radiodurans wild type and M.DraR1 mutant strains, suggesting that there might be a relevant cognate REase limiting the foreign DNA (containing "CCGCGG" motifs) uptake in D. radiodurans. A putative GIY-YIG endonuclease encoded by ORF15945P was also annotated in the genome of D. radiodurans, which may recognize and cleave the same motif based on the high similarity of its amino acid sequence with that of R.Eco29kI.
Cells lacking 4mC methylation displayed decreased genomic stability in D. radiodurans. Previous studies proposed that DNA methylation might affect the chromosome topology whereby induces conformational changes to a bacterial genome, rendering previously inaccessible genes accessible to the transcriptional machinery (Polaczek et al., 1998;Ngo et al., 2016;Beaulaurier et al., 2019). Interestingly, almost all nucleoid-associated  proteins, e.g., Dps-1, HU, UDG, and type IB topoisomerase, were found to be significantly downregulated in the M.DraR1 strain, suggesting the presence of defects in DNA compaction in the mutant strain (Blasius et al., 2008;Nguyen et al., 2009;Reon et al., 2012;Wang et al., 2012). DNA gyrase, containing GyrA and GyrB, which are involved in DNA relaxation, decatenation and supercoiling (Kota et al., 2016), was upregulated in M.DraR1, implying the role of 4mC methylation in the regulation of DNA condensation. The "CCGCGG" sites also found to be uniformly distributed across the D. radiodurans genome, indicating that they may also affect DNA compaction through epigenetic modulation of gene expression in D. radiodurans.

CONCLUSION
We characterized N 4 -cytosine DNA methylation in the extremophile D. radiodurans R1 and addressed its possible association with transcriptional expression patterns. Evidence from our study suggested that M.DraR1 is a novel α-class N 4 -cytosine MTase that specifically recognizes the target sequence "CCGCGG" and methylates the second cytosine at the N 4 position, preventing it from being digested by cognate endonucleases both in vitro and in vivo. Moreover, cells lacking 4mC DNA modifications exhibited higher spontaneous mutation frequencies and enhanced recombination and transformation efficiencies. And a large number of DEGs were observed in conjunction with the loss of the MTase enzyme. Our study has provided new insights into the properties and roles of N 4 -cytosine DNA methylation in the maintenance of genomic integrity in D. radiodurans.

Strains and Growth Conditions
All strains and plasmids used in this study are listed in Supplementary Table S1. E. coli strains were grown in Luria-Bertani (LB) liquid media or on agar (1.5% Bacto-agar) plates supplemented with appropriate antibiotics at 37 • C with aeration. All D. radiodurans strains were grown in tryptone glucose yeast extract (TGY) liquid media or on agar plates (0.5% tryptone, 0.1% glucose, and 0.3% yeast extract) supplemented with appropriate antibiotics at 30 • C with aeration.

SMRT-Seq
Genomic DNA was extracted using a QIAamp DNA mini kit (Qiagen, Valencia, CA, United States) following the manufacturer's protocol. The DNA quality was determined using gel electrophoresis and a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States). SMRT sequencing and de novo assembly were carried out as described previously (Hua and Hua, 2016). DNA modification detection and motif analysis were performed using the standard settings in the "RS Modification and Motif Analysis.1" protocol included in SMRT Portal version 2.2.0. 1 Modified sites were then identified through kinetic analysis of the aligned DNA sequence data and grouped into motifs using Motif Finder (v1)2 Blow et al., 2016). The SEQWARE computer resource was used to annotate R-M system genes (Blow et al., 2016). The results were stored into REBASE (Roberts et al., 2015). PacBio motif files and sequencing data were submitted to the NCBI database (Table 1).

M.DraR1 Mutant Strain Construction and Complementation
Tripartite ligation and double-crossover recombination methods were used for gene mutation as described previously with slight modifications (Dai et al., 2018). Briefly, upstream and 1 http://www.pacb.com/devnet/code.html downstream fragments of the target gene were PCR-amplified, digested and ligated to a streptomycin resistance cassette. The tripartite ligated fragment were PCR-amplified using P1/P4 primers and transformed into D. radiodurans competent cells. Mutant colonies were screened on TGY plates containing 10 µg/mL streptomycin (TS plate) and confirmed by PCR analysis and DNA Sequencing. The mutant strain was named M.DraR1. For gene complementation, the full-length M.DraR1 gene fragment was PCR-amplified using M.DraR1-F and M.DraR1-R primers and ligated into the shuttle vector pRADK to construct the complementation plasmid pRAD-M.DraR1. After confirmation by sequencing using primers pRAD-F and pRAD-R, the plasmid was subsequently transformed into M.DraR1 complement cells and selected on TGY plates containing streptomycin and 4 µg/mL chloramphenicol (TSC plate). The resulting strain was designated M.DraR1/pk-M.DraR1. All primers are listed in Supplementary Table S2.

Quantification of Modified Bases in gDNA by UHPLC-QQQ-MS/MS
The gDNA (2-4 µg) in 26 µL of nuclease-free H 2 O was denatured at 100 • C for 5 min, chilled on ice for 2 min, and digested by 1 µL nuclease P1 (1U/µL, Wako United States) in 10 mM NH 4 OAc (adding 3 µL 100 mM NH 4 OAc, pH 5.3) at 42 • C overnight. It was followed by addition of 3.4 µL NH 4 HCO 3 (1 M) and 1 µL of phosphodiesterase I from Crotalus adamanteus venom (0.001 U, Sigma-Aldrich) at 37 • C for 2 h and finally by addition of 1 U of alkaline phosphatase from E. coli (Takara) at 37 • C for 2 h. Digested DNA was diluted twofold with nucleasefree H 2 O. The diluted solution was filtered through a 0.22 µm filter and 5 µL was injected for liquid chromatography with tandem mass spectrometry (Greer et al., 2015;Liu et al., 2016).
2 -Deoxyguanosine (dG, D7145), Thymidine (dT, T9250), 2 -Deoxyadenosine (dA, D7400), 2 -Deoxycytidine (dC, D3897), and N 6 -Methyl-2 -deoxyadenosine (m 6 dA or 6mA, M2389) were purchased from Sigma. 2 -Deoxy-5-methylcytidine (m 5 dC or 5mC, D3610) was purchased from Tokyo Chemical Industry (TCI, Japan). 2 -Deoxy-N 4 -methylcytidine (m 4 dC or 4mC) was a generous gift from Prof. Jianzhao Liu (Zhejiang University, China). Nucleosides were separated by reversed-phase UHPLC on a C18 column with online mass spectrometry detection using a Waters TQ MS triple-quadrupole LC mass spectrometer in positive electrospray ionization mode (Waters, Milford, MA, United States). The elution program was shown as following: the C18 column was washed using ddH 2 O and equilibrated with 98% solvent A (1 Formic acid/ddH 2 O) and 2% solvent B (1 Formic acid/methanol); the content of B was increased gradually from 2 to 100% in 5 min after injection, and subsequently continue to elute 2.5 min using 100% B; and then the column was treated circularly as the first two steps did before the next injection. Nucleosides were quantified using nucleoside precursor ion to base ion mass transitions of 228.2-112.1 for dC, 242.2-126.1 to separate m 4 dC and m 5 dC by retention time, 266.1-150.0 for m 6 dA and 252.1/136.0 for dA (Yu et al., 2015;Liu et al., 2016). Ratio quantification was performed using calibration curves obtained from nucleoside standards run at the same time.

Methyltransferase Cloning
The M.DraR1 gene was amplified by PCR using the primers M.DraR1-pRRS-F and M.DraR1-pRRS-R and cloned into the pRRS plasmid as described previously . The pRRS-M.DraR1 plasmid was transformed into E. coli DH5α cells and cultured on ampicillin-resistant LB plates for 20 h at 37 • C. Clones were selected and confirmed by DNA sequencing using the sequencing primers pRRS-F and pRRS-R.
The constructed pRRS-M.DraR1 vector was transferred to E. coli ER2796 cells. Plasmid DNA was isolated and the methylation status was assessed by digestion with HindIII and a relevant cognate REase SacII. The unmethylated control substrate was produced by PCR amplification using primers pRRS-F and pRRS-F-R. Null pRRS was analyzed as well. Additional REases, including PvuII, HhaI, HpaII, and HaeIII were also used to digest the substrates. Samples were separated by 5% native polyacrylamide gel electrophoresis (PAGE). Gels were stained using ethidium bromide (EB) and bands were visualized using an Automatic Digital Gel Image Analysis System (Tannon 3500, Shanghai, China). The pRRS plasmid and E. coli ER2796 host strain were gifts from Sir Richard J. Roberts (New England Biolabs, Ipswich, MA, United States). All primer sequences are listed in Supplementary Table S2. All restriction enzymes were purchased from New England Biolabs, United States.

Expression and Purification of M.DraR1 MTase
The M.DraR1 gene was amplified and cloned into a modified pET28a expression vector (Zhao et al., 2015). After confirmed by DNA sequencing, the constructed plasmid was transformed into ER2566 strain. Protein expression was induced at 16 • C by adding isopropyl-β-D-thiogalactopyranoside (IPTG) to a final concentration of 0.4 mM.
The M.DraR1 protein was purified using GE Akta Pure 25 system (GE Healthcare, United States). Briefly, cells resuspended in lysis buffer (20 mM Tris-HCl, pH 8.0; 500 mM NaCl; 5% (w/v) glycerol; 3 mM β-mercaptoethanol and 9 mM imidazole) were lysed by sonication and clarified by centrifugation at 20,000 g for 30 min at 4 • C. The supernatant was purified using a HisTrap HP column (GE Healthcare) equilibrated with buffer A (20 mM Tris-HCl pH 8.0; 500 mM NaCl; 5% (w/v) glycerol and 3 mM β-mercaptoethanol), washed with 40 mM imidazole and eluted with 250 mM imidazole. After TEV tag removal using TEV protease, the protein was loaded onto an MBPTrap HP column (GE Healthcare) to remove His-MBP and uncleaved proteins. The flow-through fractions were collected, desalted using a HiPrep TM 26/10 column (GE Healthcare) with buffer B (20 mM Tris-HCl, pH 8.0; 100 mM NaCl; 1 mM dithiothreitol; 5% (w/v) glycerol and 1mM EDTA) and loaded onto a Heparin HP column (GE Healthcare) pre-equilibrated with the same buffer. Fractions containing M.DraR1 protein were eluted using a linear NaCl gradient from 100 to 1000 mM. Finally, the protein was purified using a Superdex 200 10/300 GL column (GE Healthcare) with buffer C (20 mM Tris-HCl, pH 8.0 and 200 mM KCl). The purified protein was assessed by sodium dodecyl sulfate-PAGE, Western blotting and matrix-assisted laser desorption/ionization time of flight mass spectrometry.

Electrophoresis Mobility Shift Assays
Electromobility shift assays were performed as previously described with some modifications (Cheng et al., 2016). Briefly, 20 µL reaction mixtures containing 100 µM double-stranded DNA (generated by annealing the oligonucleotides listed in Supplementary Table S2) and various concentrations of purified M.DraR1 (0, 100, 200, 300, 400, 500, 600, and 1000 µM) in binding buffer (20 mM Tris-HCl, pH 8.0 and 100 mM KCl) with excess or no SAM, were incubated for 60 min at 4 • C. Samples were analyzed by 12% native PAGE gel electrophoresis and bands were visualized by EB staining. All DNA oligonucleotides were purchased from Sangon.

In vitro M.DraR1 Protection Assays
For in vitro methylation, 1 µg unmethylated pRRS-M.DraR1 plasmid, 2 µg DNA containing three "CCGCGG" sites (PCR products of λDNA and gDNA of D. radiodurans, 2,836 and 1,318 bp, respectively, primers listed in Supplementary Table S2) and λDNA were incubated with 1 µM purified M.DraR1 at 30 • C for 1 h in reaction buffer (20 mM Tris-HCl, pH 8.0; 100 mM KCl; 0.1 mM EDTA; 3 mM β-mercaptoethanol and 80 µM SAM). Methylated products were then purified and digested as described above. Finally, the digested products were analyzed by 1% agarose or 5% native PAGE gel electrophoresis as described above.

Growth Curve Test
Growth curve and phenotypic assays were performed as previously described (Cheng et al., 2015). D. radiodurans and derivative strains were grown to early exponential phase (OD 600 ≈ 1.0) and chilled on ice. To assess growth rate, 0.5 mL aliquots were resuspended in 50 mL fresh TGY media and incubated at 30 • C with shaking at 220 rpm. The OD 600 was measured every 4 h for a total incubation period of 40 h. Growth curves were created from at least three independent measurements using GraphPad Prism (GraphPad Software, Inc., La Jolla, CA, United States).

Transformation, Recombination, and Spontaneous Mutation Frequency
A total of 0.2 µg of each of the shuttle vector pRADK, non-methylated and methylated pRADKm (containing one "CCGCGG" site inserted by site-directed mutagenesis using primers pRADKm-F and pRADKm-R and confirmed by sequencing using primer pRADKm-seq) was transferred into D. radiodurans strains. After incubation, cultures were plated onto TGY and TGY containing chloramphenicol (TC) plates. The transformation frequency was determined by dividing the number of transformant colonies growing on TC plates by the total number of viable cells on TGY plates (Cheng et al., 2015). We also amplified a donor DNA fragment containing a kanamycin resistance gene from crtB mutant (Han et al., 2017). A total of 0.5 µg this fragment was used to measure the recombination capacity of the M.DraR1 mutant. Approximately 10 9 freshly cultured cells were harvested at the early exponential phase, resuspended and plated on TGY plates containing rifampicin to determine the number of mutants and on native TGY plates to determine the total number of viable cells (Cheng et al., 2015). Three replicates were performed for each strain.

RNA Isolation and Sequencing (RNA-Seq)
Total RNA of two biological replicates (OD 600 = 0.5) from each group was isolated using the TRIzol method according to the manufacturer's protocol (Ambion, Foster City, CA, United States). RNA degradation and contamination were assessed using 1% agarose gels. RNA concentration and RNA integrity were assessed using an RNA Nano 6000 Assay kit and Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). A total of 3 µg RNA per sample was used as the input material for library preparation. Index codes were added to attribute sequences to each sample. The library preparations were sequenced on an Illumina Hiseq platform and paired-end reads were generated. The raw reads were deposited in the Sequence Read Archive (SRA) database of NCBI (PRJNA532990).

RNA-Seq Data Analysis
The index of the reference genome (downloaded directly from the NCBI) was built and clean reads were aligned to the reference genome by Bowtie2-2.2.3 (Langmead and Salzberg, 2012). HTSeq v0.6.1 was used to count the number of reads mapped to each gene. Differential expression analysis was performed using the DEGs analytical method (1.18.0), which provides statistical routines for determining differential expression in digital gene expression data using a negative binomial distribution model (Anders and Huber, 2010). GO enrichment analysis of DEGs was performed using the GOseq R package, which corrected for gene length bias. GO terms with corrected p-values < 0.05 were considered significantly enriched (Young et al., 2010). We also test the enrichment of DEGs in Kyoto Encyclopedia of Genes and Genomes pathways (Mao et al., 2005).

Experimental Validation
The reliability of the RNA-seq data was verified using 13 DEGs with different fold-change values (seven up-regulated and six down-regulated genes). Total RNA was extracted and cDNA was amplified using SYBR Premix Ex Taq TM (TaKaRa). RT-qPCR assays were performed using the Mx3005 TM P Real-time detection system (Agilent, United States) and DR_1343 was used to normalize expression. RT-qPCR primers are listed in Supplementary Table S2.

DATA AVAILABILITY
The RNA−seq data from this study has been submitted to the NCBI Gene Expression Omnibus with an accession number of GSE134934 (https://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE134934).

AUTHOR CONTRIBUTIONS
YH conceived the project. SL, YZ, and BT designed the experiments and wrote the manuscript. SL, HL, and JC carried out the construction, purification, mass spectrometric analysis, and biochemical experiments. SD, JH, and SM carried out the mutants and phenotypic assays. SL, LW, and XH carried out the SMRT-seq, RNA-seq, and data analysis. HX and SL performed the enzymatic activity assays. All authors took part in the data analysis and approved the final version of the manuscript to be published.