text:
Characterizing the DNA Methyltransferases of Haloferax volcanii via Bioinformatics, Gene Deletion, and SMRT Sequencing

DNA methyltransferases (MTases), which catalyze the methylation of adenine and cytosine bases in DNA, can occur in bacteria and archaea alongside cognate restriction endonucleases (REases) in restriction-modification (RM) systems or independently as orphan MTases. Although DNA methylation and MTases have been well-characterized in bacteria, research into archaeal MTases has been limited. A previous study examined the genomic DNA methylation patterns (methylome) of the halophilic archaeon Haloferax volcanii, a model archaeal system which can be easily manipulated in laboratory settings, via single-molecule real-time (SMRT) sequencing and deletion of a putative MTase gene (HVO_A0006). In this follow-up study, we deleted other putative MTase genes in H. volcanii and sequenced the methylomes of the resulting deletion mutants via SMRT sequencing to characterize the genes responsible for DNA methylation. The results indicate that deletion of putative RM genes HVO_0794, HVO_A0006, and HVO_A0237 in a single strain abolished methylation of the sole cytosine motif in the genome (Cm4TAG). Amino acid alignments demonstrated that HVO_0794 shares homology with characterized cytosine CTAG MTases in other organisms, indicating that this MTase is responsible for Cm4TAG methylation in H. volcanii. The CTAG motif has high density at only one of the origins of replication, and there is no relative increase in CTAG motif frequency in the genome of H. volcanii, indicating that CTAG methylation might not have effectively taken over the role of regulating DNA replication and mismatch repair in the organism as previously predicted. Deletion of the putative Type I RM operon rmeRMS (HVO_2269-2271) resulted in abolished methylation of the adenine motif in the genome (GCAm6BN6VTGC). Alignments of the MTase (HVO_2270) and site specificity subunit (HVO_2271) demonstrate homology with other characterized Type I MTases and site specificity subunits, indicating that the rmeRMS operon is responsible for adenine methylation in H. volcanii. Together with HVO_0794, these genes appear to be responsible for all detected methylation in H. volcanii, even though other putative MTases (HVO_C0040, HVO_A0079) share homology with characterized MTases in other organisms. We also report the construction of a multi-RM deletion mutant (ΔRM), with multiple RM genes deleted and with no methylation detected via SMRT sequencing, which we anticipate will be useful for future studies on DNA methylation in H. volcanii.

Introduction
In bacteria and archaea, DNA methylation by DNA methyltransferases (MTases) has many roles. MTases are commonly associated with restriction-modification (RM) systems, in which the MTase functions alongside a cognate restriction endonuclease (REase). The REase will target the same sites of DNA as the MTase and cleave those that are unmethylated, whereas methylated motifs will be disregarded. RM systems function in self-recognition, allowing the host to differentiate between its own methylated DNA and potentially harmful foreign unmethylated DNA, which can then be recognized and digested by the REase [1,2]. RM systems have also been characterized as Toxin/Antitoxin systems that lead to addiction through post-segregational killing in case the methylation activity decays in a cell from which the RM system has been lost so that unmethylated restriction sites become exposed to the remaining restriction enzyme activity [3,4].
DNA methylation occurs at adenine or cytosine bases, resulting in one of three possible types of methylation: N6-methyladenine (6mA), N4-methylcytosine (4mC), or C5-methylcytosine (5mC) [5]. Methylation is catalyzed by DNA MTases, which interact with the cofactor S-adenosyl methionine (AdoMet) to transfer a methyl group to a nucleotide base of a DNA molecule. MTases typically consist of three major domains: an AdoMet binding domain which interacts with AdoMet to obtain the methyl group, a target recognition domain (TRD) which recognizes a short sequence of DNA to be targeted for methylation, and a catalytic domain which transfers the methyl group from AdoMet to the targeted nucleotide [6]. MTases can be categorized based on the type of methylation they perform (6mA, 4mC, and 5mC), with the 6mA and 4mC MTases being more similar to each other than to 5mC MTases [7,8].
MTases have been further categorized into subtypes based on the order of several conserved motifs that make up the primary domains of the MTases. The 6mA and 4mC MTases can be classified into six categories (α, β, γ, δ, ε, and ζ) based on the N-terminal to C-terminal order of conserved motifs X, I-III (AdoMet binding motifs), IV-VIII (catalytic motifs), and the TRD [8,9]. The occurrence of signature AdoMet binding motif DPPY and catalytic motif FXGXG (abbreviated FGG) can also be used to categorize these MTases [10]. The 5mC MTases have a different set of motifs that can be used to identify them [11,12].
There are four major types of characterized RM systems [13,14]. Type I RM systems consist of pentamer complexes with two REase (R) subunits, two MTase (M) subunits, and one site specificity (S) subunit containing two tandem TRDs which recognize bipartite target sites. When the complex comes across a target site, it will either methylate the site if it is methylated on one strand (hemimethylated) or will cleave the DNA several bases upstream or downstream from the site if it is unmethylated on both strands [15,16]. Type II RM systems include MTases and REases which operate independently and target the same sites of DNA [17]. Many subgroups of Type II RM systems have been categorized, such as Type IIG which consists of independent RM enzymes capable of both MTase and REase activity [18]. In Type III RM systems, a REase subunit (Res) and MTase subunit (Mod) work together in a two-component complex, with the Mod subunit containing the TRD which recognizes the target site of the system [19]. Type IV RM systems consist only of REases, and these REases cleave methylated target sites instead of unmethylated sites [20].
MTases can also occur independently in bacteria and archaea without cognate REases. These MTases, known as orphan MTases, typically provide important functions for their host organisms [21]. In Escherichia coli, for example, the orphan adenine MTase Dam is involved in coordinating timing of DNA replication by methylating GATC sites at the origin of replication which are also targets of binding for SeqA in a hemimethylated state [22]. When SeqA binds to hemimethylated GATC sites of the origin after replication has occurred, DnaA is prevented from binding to the origin and re-initiating DNA replication [23,24]. Dam methylation is also important in the methyl-directed mismatch repair (MMR) system in E. coli, in which the complex binds to a closely located methylated GATC site on the old strand in order to target and cleave the mismatched base on the new strand [25][26][27]. In Caulobacter crescentus, the orphan adenine MTase CcrM is involved in regulating the expression of genes like ctrA, which are essential for cell cycle regulation [28,29]. Orphan MTases can also protect the host from parasitic RM systems by mimicking the methylation of the invader, such as orphan cytosine MTase Dcm in E. coli which methylates the same sites as RM system EcoRII and prevents degradation of the genome by the invading system [30]. Orphan MTases are more common among the bacteria and more well-conserved within a genus than RM-associated MTases, likely due to orphan MTases performing important roles within their hosts [31,32]. DNA methylation has been well-studied in bacterial organisms. However, research has not been as extensive in the archaea, which have focused on characterizing methylation and a few RM systems primarily in thermophilic organisms [33][34][35][36]. A previous study [37] examined the genomic methylation patterns (methylome) of the halophilic archaeal organism Haloferax volcanii, a member of Class Halobacteria which are often referred to as haloarchaea, as a model for examining DNA methylation and RM systems in archaea, due to its well-established genetic system which allows it to be easily manipulated in lab settings [38,39]. In the study, the methylome of H. volcanii was sequenced via single-molecule real-time (SMRT) sequencing developed by Pacific Biosciences (PacBio) [40]. H. volcanii was observed to have two types of motifs methylated throughout its genome: C m4 TAG and GCA m6 BN 6 VTGC. The study also demonstrated that deletion of one of the putative RM genes (HVO_A0006) resulted in an alteration in the adenine motif, which was surprising considering that the gene is not a MTase but instead encodes an REase family protein [37]. In this follow-up study, we aim to characterize the MTases of H. volcanii through bioinformatics and gene deletions of the various predicted RM genes in the genome and sequence the methylomes of the deletion mutants via SMRT sequencing. We will also describe the production of an RM null mutant without a methylated genome, which we anticipate will be useful in future research of DNA methylation in the archaea.

Strains and Growth Conditions
All strains and plasmids used in this study are listed and described in Table 1. Strains of H. volcanii were grown at 42 • C while shaking at 200 rpm using either rich medium (Hv-YPC) or selective rich medium (Hv-Ca) developed by Allers et al. [41] and outlined in the Halohandbook [42]. For ∆pyrE2 strains, media was supplemented with uracil (50 µg/mL) and 5-fluoroorotic acid (50 µg/mL) as needed. Strains of E. coli were grown at 37 • C while shaking at 200 rpm in either Lysogeny Broth (LB) or S.O.C. medium (Clontech, Mountain View, CA, USA). Ampicillin (100 µg/mL) and X-gal (20 µg/mL) were added to the media when needed.

Deletion of Annotated Restriction Modification Genes
Putative RM genes in H. volcanii were identified from New England BioLabs Restriction Enzyme Database (REBASE) [10] and National Center for Biotechnology Information (NCBI) ( Table 2). These genes were deleted in H. volcanii strain H1206 utilizing a method developed by Blaby et al. [39] that uses the In-Fusion HD Cloning Kit (Clontech). Primers were designed to construct deletion plasmids of putative RM genes and are listed in Table 3. These deletion plasmids were then used to transform H. volcanii H1206 and its derivatives via the polyethylene glycol (PEG)-mediated transformation protocol outlined in the Halohandbook [42]. Transformed cell cultures were plated on Hv-Ca and incubated at 42 • C for 5-7 days. Pop-ins were detected via a colony PCR screen using external deletion primers and visualized via gel electrophoresis. Confirmed pop-ins were then plated on Hv-Ca with 50 µg/mL 5-fluoroorotic acid (5-FOA) and 50 µg/mL uracil to pop-out genes of interest. Successful pop-outs were identified via PCR screen as performed for detecting pop-ins. Final deletion strains obtained though this method are listed in Table 1.

DNA Purification for Single-Molecule Real-Time Sequencing
In order to extract DNA from the H. volcanii deletion mutants for SMRT sequencing, 40 mL of cell cultures in late log to early stationary phase (optical density (OD 600 ) =~0.8-1) were pelleted and lysed by resuspension in 10 mM Tris-HCl buffer (pH 8.0). The lysates were then treated with proteinase K (50 µg/mL final concentration) and incubated overnight at 37 • C to hydrolyze the proteins, after which the DNA was extracted via ethanol precipitation. Performing three rounds of phenol-chloroform and two rounds of chloroform extractions purified the DNA further. A final ethanol precipitation was then performed on the remaining DNA, and the samples were purified of RNA via Agencourt AMPure XP beads (Beckman Coulter, Brea, CA, USA). The 260/280 ratio, 260/230 ratio, and DNA concentration of each sample was quantified via Nanodrop and Qubit dsDNA BR assay (Invitrogen, Eugene, OR, USA).

Single-Molecule Real-Time Sequencing
The DNA samples extracted from the H. volcanii deletion mutants were analyzed via PacBio SMRT sequencing in order to determine the methylomes of the strains. The samples were submitted to the Keck Sequencing Facility of the Yale School of Medicine for SMRT sequencing analysis. A detailed outline of the SMRT sequencing strategy can be found in the PacBio manual "Detecting DNA Base Modifications: SMRT Analysis of Microbial Methylomes" [46]. Libraries of 0.25 to 3 kb were constructed for each strain using an estimated input size of 4 Mb, and were each sequenced in one SMRT cell, resulting in coverage of~150x for ∆RM,~400x for ∆HVO_0794 ∆HVO_A0006 ∆HVO_A0237, and~120x for ∆rmeRMS. The SMRT Portal program "RS_Modification_and_Motif_Analysis.1" was used under default settings to determine the modified bases and motifs in ∆RM. The modified bases and motifs in the other strains were identified using the same SMRT Portal program, but with the ∆RM analysis results used as a control. All analyses used the H. volcanii DS2 genome as the reference sequence [47].

Bioinformatics Analysis
Homologs of the putative RM proteins in H. volcanii DS2 were discovered via protein BLAST (blastp) [48] and position-specific iterative BLAST (PSI-BLAST) of the non-redundant protein database on NCBI as well as translated nucleotide BLAST (tblastn) of the NCBI Halobacteria genome database (taxid 183963) (E-value cutoff of 1e −4 ). Homologs were also identified using the REBASE database of RM genes [10]. Alignments of identified homologs were performed using Clustal X2 [49]. Protein domain architecture and sequence features, including identification of Structural Classification of Proteins (SCOP) superfamilies, were analyzed using InterProScan [50].
Homologs to the CTAG modification methyltransferase in H. volcanii DS2 (ADE02643) in completely sequenced halobacterial genomes were identified using the NCBI's blast site for microbial genomes, selecting completely sequenced halobacterial genomes and the tblastn search algorithm. The list of completely sequenced genomes did not completely correspond to the genomes searched through the NCBI's web interface; therefore, the absence of a homolog in a genome was confirmed through a targeted tblastn search. The one additional homologous gene identified in this step was added to the phylogenetic analysis. Matching nucleotide sequences were retrieved, translated into protein and aligned using muscle [51] as implemented in Seaview [52], and used for phylogenetic reconstruction using PhyML [53] with the following parameters: LG substitution model, 100 bootstrap samples, 4 substitution rate categories with estimated Gamma distribution parameter, and estimated fraction of invariant sites, and a tree topology search using both Nearest Neighbor Interchange and Subtree Pruning and Regrafting.
CTAG and GATC frequency and cumulative occurrence of these motifs were calculated with an in house Perl (Practical Extraction and Report Language) script.

Haloferax volcanii Growth Experiments
Haloferax volcanii strains H26 and ∆RM (Table 1) were grown in Hv-YPC medium to mid-log phase (OD 600 =~0. 6-0.8). The cell cultures were then diluted in Hv-YPC to an OD 600 of~0.01 and were each distributed into 24 wells of a 96-well plate, with each well receiving 200 µL of culture. One well on the plate received 200 µL of Hv-YPC to be used as a blank reading. The 96 well plate was then covered with sealing tape and inserted into a Multiscan FC plate reader (Fisher Scientific, Waltham, MA, USA), which recorded the OD 620 of each well every hour for 72 h while incubating and shaking the plate at 42 • C.

Bioinformatics Analysis Supports Identification of HVO_0794 as a Chromosomal 4mC CTAG MTase
The putative 4mC CTAG MTase HVO_0794 was analyzed bioinformatically. A blastp analysis identified a homolog to the enzyme in Methanothermobacter thermautotrophicus named M.MthZI (GenBank CAA48447) which has been experimentally characterized as a 4mC CTAG MTase [33]. Two other homologs were also identified via blastp that were also experimentally characterized via unpublished work according to REBASE: M.BfaI (GenBank ADQ20483) in Bacteroides fragilis and M.MjaI (GenBank AAB98988) in Methanocaldococcus jannaschii DSM 2661. These homologs are similar in size to HVO_0794, ranging between 303 to 364 amino acids in length. Also, these enzymes are classified as Type II, subtype β 4mC MTases on REBASE, as is HVO_0794. A multiple sequence alignment ( Figure 1) of HVO_0794 with these homologs and homolog M.HsaR1I (GenBank CAP14114) from Halobacterium salinarium R1 indicate significant sequence similarity is shared in the N-terminal and central regions of the amino acid sequences. This region of sequence similarity belongs to the S-adenosyl-L-methionine-dependent methyltransferase superfamily domain SSF53335 identified by InterProScan in the SCOP database. Signature N4-methyltransferase motifs PR00508 from the protein motif database PRINTS were also identified in the region via InterProScan (data not shown). A closer examination of the alignment revealed the presence of motifs I-X identified in M.MthZI and other 4mC MTases by Bujnicki and Radlinska [9]. These motifs are present in the alignment in the order of N-III-IV-V-VI-VII-VIII-VIII'-IX-X-I-II-C which is indicative of subtype β MTases [8]. The signature AdoMet binding motif DPPY and catalytic motif FGG are also fully conserved in the alignment (DPPY conserved here as SPPY). The FGG motif also occurs before the DPPY motif in the alignment, a motif order observed in subtype β MTases according to REBASE [10]. Overall, these results support the identification of HVO_0794 as a 4mC CTAG MTase of the subtype.
In a search of all completely sequenced halobacterial genomes available on 12 June 2017, homologs that group with HVO_0794 with high statistical support were identified in 37 (88%) of the completely sequenced genomes. Genomes with the homolog present are included in Figure 2. Homologs that grouped with the H. volcanii CTAG enzyme with high support were absent in Halorubrum lacusprofundi ATCC 49239, Halorubrum trapanicum, Haloquadratum walsbyi C23, Haloquadratum walsbyi DSM 16790, and Halopenitus persicus.
The GATC motif is methylated in many organisms and this methylation was shown to play a role in regulating the start of replication and in mismatch repair of newly synthesized DNA strands [23][24][25][26][27]. CTAG and GATC motifs occur throughout the H. volcanii genome. In contrast to other Halobacteria, both motifs show localized areas of higher concentrations within the genome. The H. volcanii chromosome possesses several origins of replication [54], one of these (oriC2) is associated with an increased concentration of CTAG and GATC motifs ( Figure 3).
In E. coli and other organisms where GATC methylation facilitates recognition of the newly synthesized DNA strand during mismatch repair, the GATC motif occurs with higher frequency as compared to the CTAG motif (22 times in E. coli, 46 times in H. trapanicum, see Table 4). In H. volcanii this ratio is only 2.8 (see Figure 2 and Table 4). Ratios below 5 were found in other Haloferax species, Haloarcula sp., Natronomonas pharaonic, Halobacterium hubeiense strain JI20-1, and Halobacterium sp. DL1; whereas Halobacterium salinarum R1 has a ratio above 20 (Table 4). In Haloferax spp. this drop in relative GATC frequency is due to a dropin frequency of the GATC ( Table 4). The CTAG motif actually occurs less frequently in H. volcanii (0.24 times per 1000 nucleotides) than in other halobacterial chromosomes (average ± standard deviation in all completely sequenced chromosomes is 0.42/kb (±0. 19).  . Identified N4-cytosine methyltransferase motifs I-X [9] are highlighted in blue (representing S-adenosyl methionine (AdoMet) binding motifs), purple (representing DNA binding motif), and yellow (representing catalytic motifs). Red boxes are used to identify the signature DPPY and FGG motifs. The SCOP superfamily domain S-adenosyl-L-methionine-dependent methyltransferase domain SSF53335 is highlighted in green throughout the alignment. Clustal X2 shading and marking of amino acids is included in the alignment.   The phylogeny was rooted using more divergent haloarchaeal and methanomicrobial homologs. Genomes with a chromosome wide GATC to CTAG ratio below five are given in red, those with a GATC to CTAG ratio between 5 and 14 are given in orange. In addition, those with a ratio above 20 are given in black. Note that only few groups are well supported, including the Haloferax and Haloarcula genera.
The GATC motif is methylated in many organisms and this methylation was shown to play a role in regulating the start of replication and in mismatch repair of newly synthesized DNA strands [23][24][25][26][27]. CTAG and GATC motifs occur throughout the H. volcanii genome. In contrast to other Halobacteria, both motifs show localized areas of higher concentrations within the genome. The H. volcanii chromosome possesses several origins of replication [54], one of these (oriC2) is associated with an increased concentration of CTAG and GATC motifs (Figure 3). The phylogeny was rooted using more divergent haloarchaeal and methanomicrobial homologs. Genomes with a chromosome wide GATC to CTAG ratio below five are given in red, those with a GATC to CTAG ratio between 5 and 14 are given in orange. In addition, those with a ratio above 20 are given in black. Note that only few groups are well supported, including the Haloferax and Haloarcula genera. In E. coli and other organisms where GATC methylation facilitates recognition of the newly synthesized DNA strand during mismatch repair, the GATC motif occurs with higher frequency as compared to the CTAG motif (22 times in E. coli, 46 times in H. trapanicum, see Table 4). In H. volcanii this ratio is only 2.8 (see Figure 2 and Table 4). Ratios below 5 were found in other Haloferax species, Haloarcula sp., Natronomonas pharaonic, Halobacterium hubeiense strain JI20-1, and Halobacterium sp. DL1; whereas Halobacterium salinarum R1 has a ratio above 20 (Table 4). In Haloferax spp. this drop in relative GATC frequency is due to a dropin frequency of the GATC ( Table 4). The CTAG motif actually occurs less frequently in H. volcanii (0.24 times per 1000 nucleotides) than in other halobacterial chromosomes (average ± standard deviation in all completely sequenced chromosomes is 0.42/kb (±0.19).

Bioinformatics Analysis Supports Identification of RmeM as a Type I 6mA MTase and RmeS as a Type I Specificity Subunit on the Chromosome
The putative Type I 6mA MTase RmeM and its cognate specificity subunit RmeS were also analyzed via bioinformatics. Tblastn of the RmeM sequence against the database of Halobacteria genomes in NCBI (taxid 183963) showed this MTase to be relatively rare in this Order, as we retrieved significant hits in 19 out of 181 (10.5%) genomes of Halobacteria, and 3 out of 42 (7.1%) of the fully sequenced genomes. Blastp analysis of RmeM indicated that it is homologous to M.EcoKI (GenBank P08957), a well-characterized Type I 6mA MTase in E. coli [55,56]. A homolog to RmeM was also identified in Bacillus cereus ATCC 10987 (M.BceSVI; GenBank AAS39772), which has been characterized via SMRT sequencing [57], as well as in Methanoregula boonei 6A8 (Mboo_1031; GenBank ABS55549) which was characterized via unpublished research according to REBASE. These enzymes are 477 to 529 amino acids in length, which is similar to RmeM in size. The classification on REBASE for these enzymes also matches RmeM (Type I, subtype γ MTase). A multiple sequence alignment (Figure 4) of these identified homologs, along with a homolog identified in Caldanaerobacter subterraneus subsp. tengcongensis MB4 (Tte_1547; GenBank AAM24756) and in Halobacterium salinarum NRC-1 (M.HspNI; GenBank AAG18733) indicated sequence similarity is shared throughout the alignment. The same SCOP superfamily domain SSF53335 present in the HVO_0794 homologs was also identified in these homologs and spans most of the alignment. InterProScan also identified PFAM protein family database domain PF12161, an N-terminal domain present in Type I MTases which affects the affinity of the MTase for hemimethylated DNA [58]. Another PFAM domain was also detected: PF02384, which is a 6mA MTase domain found in Type I MTase enzymes. The catalytic signature motif FGG is conserved (AGG in this alignment), although the third residue of the motif is not well-conserved in M.BceSVI and M.EcoKI. The AdoMet binding signature motif DPPY is also conserved in the alignment as NPP(Y/F). Both of these signature motifs are poorly conserved in M.HspNI. The order of these signature motifs in the alignment, with FGG occurring before DPPY, is typical of subtype γ MTases according to REBASE [10]. These data indicate that RmeM is a Type I, subtype γ MTase. domain was also detected: PF02384, which is a 6mA MTase domain found in Type I MTase enzymes. The catalytic signature motif FGG is conserved (AGG in this alignment), although the third residue of the motif is not well-conserved in M.BceSVI and M.EcoKI. The AdoMet binding signature motif DPPY is also conserved in the alignment as NPP(Y/F). Both of these signature motifs are poorly conserved in M.HspNI. The order of these signature motifs in the alignment, with FGG occurring before DPPY, is typical of subtype γ MTases according to REBASE [10]. These data indicate that RmeM is a Type I, subtype γ MTase.  Tblastn of the RmeS sequence against the NCBI database of Halobacteria genomes (taxid 183963) resulted in significant hits in 26 out of 181 (14.4%) genomes of Halobacteria, and 7 out of 42 (16.7%) of the fully sequenced genomes. Further analysis of RmeS via blastp and PSI-BLAST resulted in the identification of a homolog in E. coli K-12 called S.EcoKI (GenBank AAC77304), a well-characterized site specificity subunit that belongs to the same Type I RM system as M.EcoKI [59]. Blastp analysis also indicated that RmeS is homologous to Tte_1545 (GenBank AAM24754), a Type I site specificity subunit in Caldanaerobacter subterraneus subsp. tengcongensis MB4 which corresponds to the same RM system as Tte_1547 and has been structurally analyzed [60]. RmeS also shares homology with S.BceSVI (GenBank AAS39773), the site specificity subunit which belongs to the same Type I RM system as M.BceSVI in Bacillus cereus ATCC 10987 [57]. RmeS was also observed to be homologous to the putative cognate site specificity subunit of M.HspNI in Halobacterium salinarum NRC-1 (S.HspNI; GenBank AAG18734) as well as the putative site specificity subunit of Mboo_1031 in Methanoregula boonei 6A8 (Mboo_1032; GenBank ABS55550). These homologs range from 398 to 476 amino acids in length, similar to RmeS which is 410 amino acids long, and were all classified as Type I specificity subunits on REBASE. A multiple sequence alignment of these enzymes ( Figure 5) did not show high sequence conservation among the homologs. However, InterProScan revealed that these homologs all shared the same SCOP superfamily DNA methylase specificity domain SSF116734. This superfamily domain was observed to occur twice in similar regions of the homologs in the alignment: one was more N-terminal in its location and the other was more C-terminal. Within these regions, the PFAM restriction endonuclease, type I, HsdS domains PF01420 were observed to occur (not shown in the alignment), which correspond to the two target recognition domains of Type I site-specificity subunits [61]. In summary, these results indicate that RmeS is the cognate Type I site-specificity subunit of RmeM.

Bioinformatics Analysis Supports Annotation of HVO_C0040 as a 5mC MTase and HVO_A0079 as a 6mA MTase
The other two putative MTase genes in H. volcanii DS2, which are located on plasmids, were also examined bioinformatically. Putative Type II 5mC MTase HVO_C0040 is located on extrachromosomal plasmid pHV1 and is flanked by an upstream IS4 family transposase (HVO_C0039). Tblastn of the HVO_C0040 sequence against the NCBI database of Halobacteria genomes (taxid 183963) resulted in significant hits in 87 out of 181 (48.1%) of halobacterial genomes, and 13 (31%) of the fully sequenced genomes. A blastp analysis of HVO_C0040 revealed that it is homologous to M.HgiDII (GenBank CAA38941) in Herpetosiphon aurantiacus, which has been experimentally characterized as a 5mC MTase recognizing GTCGAC [62]. HVO_C0040 was also observed to share homology with M.BbrUII (GenBank ABE95799), which has been characterized in Bifidobacterium breve UCC2003 as a 5mC MTase [63]. Two other homologs identified via blastp include the putative 5mC MTase in halobacterial species Halorhabdus tiamatea SARL4B (M.Hti4BORF752P; GenBank CCQ33914). Putative 5mC MTase in Acinetobacter baumannii MAR002 (M.AbaMAR002ORF10745P; GenBank KGF60346) was also identified as a homolog via blastp. These sequences range from 347 to 415 amino acids in length, which is similar to the 406-amino acid length of HVO_C0040. These sequences are also all annotated as Type II 5mC MTases on REBASE. An amino acid alignment of these homologs ( Figure 6) indicates that sequence similarity is shared in many

Bioinformatics Analysis Supports Annotation of HVO_C0040 as a 5mC MTase and HVO_A0079 as a 6mA MTase
The other two putative MTase genes in H. volcanii DS2, which are located on plasmids, were also examined bioinformatically. Putative Type II 5mC MTase HVO_C0040 is located on extrachromosomal plasmid pHV1 and is flanked by an upstream IS4 family transposase (HVO_C0039). Tblastn of the HVO_C0040 sequence against the NCBI database of Halobacteria genomes (taxid 183963) resulted in significant hits in 87 out of 181 (48.1%) of halobacterial genomes, and 13 (31%) of the fully sequenced genomes. A blastp analysis of HVO_C0040 revealed that it is homologous to M.HgiDII (GenBank CAA38941) in Herpetosiphon aurantiacus, which has been experimentally characterized as a 5mC MTase recognizing GTCGAC [62]. HVO_C0040 was also observed to share homology with M.BbrUII (GenBank ABE95799), which has been characterized in Bifidobacterium breve UCC2003 as a 5mC MTase [63]. Two other homologs identified via blastp include the putative 5mC MTase in halobacterial species Halorhabdus tiamatea SARL4B (M.Hti4BORF752P; GenBank CCQ33914). Putative 5mC MTase in Acinetobacter baumannii MAR002 (M.AbaMAR002ORF10745P; GenBank KGF60346) was also identified as a homolog via blastp. These sequences range from 347 to 415 amino acids in length, which is similar to the 406-amino acid length of HVO_C0040. These sequences are also all annotated as Type II 5mC MTases on REBASE. An amino acid alignment of these homologs ( Figure 6) indicates that sequence similarity is shared in many regions of the amino acid sequences. Many of these regions where significant sequence similarity is observed are identified as signature 5mC MTase motifs. Three regions of sequence similarity, for example, are identified by InterProScan as PRINTS cytosine-specific DNA MTase signature domains PR00105. Other regions match the 5mC conserved signature motifs identified by Posfai et al. [11], such as FGG, PC, ENV, QRR, and YGN (conserved here as (R/L)GN). Each sequence also belongs to the SCOP S-adenosyl-L-methionine-dependent MTase superfamily domain SSF53335 identified by InterProScan. Overall, these data support the annotation of HVO_C0040 as a 5mC MTase.  Putative Type IIG 6mA MTase HVO_A0079 is located on extrachromosomal plasmid pHV4 and is flanked by a downstream IS4 family transposase (HVO_A0080). Tblastn of the HVO_A0079 sequence against the NCBI database of halobacterial genomes (taxid 183963) resulted in significant hits in 101 out of 181 (59.1%) genomes of Halobacteria, and 21 (50%) of the fully sequenced genomes, indicating a higher prevalence in the Order than the RmeM or RmeS homologs, but not as high as HVO_0794 homologs. Blastp analysis of putative Type IIG 6mA MTase HVO_A0079 identified a number of homologs to the protein, including RM.Aco12261II (GenBank ADE57453), a Type IIG 6mA MTase in Aminobacterium colombiense DSM 12261 which has been identified via SMRT sequencing as targeting the motif CCRGA m6 G [32]. HVO_A0079 was also observed to share homology with RM.Fla104114II (GenBank BAV07385), a Type IIG 6mA MTase characterized via unpublished SMRT sequencing data on REBASE. Blastp also determined that HVO_A0079 shared homology with putative Type IIG 6mA MTases in Halorubrum californiensis DSM19288 (C463_0072; GenBank ELZ48543) and in Halophilic archaeon DL31 (RM.HarDL31ORF105P; GenBank AEN07377). These proteins range between 1117 to 1185 amino acids in length, which is similar to the 1088 amino acid length of HVO_A0079. These homologs are also all annotated on REBASE as Type IIG 6mA subtype α RM proteins, with the exception of C463_0072, which is not present in REBASE. A multiple GenBank CAA38941). The protein motif database PRINTS cytosine-specific DNA methyltransferase signature domains PR00105 are highlighted in yellow. Red boxes are used to identify signature FGG, PC, ENV, QRR, and YGN motifs described in Posfai et al. [11]. The SCOP superfamily domain S-adenosyl-L-methionine-dependent methyltransferase domain SSF53335 is highlighted in green throughout the alignment. Clustal X2 shading and marking of amino acids is included in the alignment.
Putative Type IIG 6mA MTase HVO_A0079 is located on extrachromosomal plasmid pHV4 and is flanked by a downstream IS4 family transposase (HVO_A0080). Tblastn of the HVO_A0079 sequence against the NCBI database of halobacterial genomes (taxid 183963) resulted in significant hits in 101 out of 181 (59.1%) genomes of Halobacteria, and 21 (50%) of the fully sequenced genomes, indicating a higher prevalence in the Order than the RmeM or RmeS homologs, but not as high as HVO_0794 homologs. Blastp analysis of putative Type IIG 6mA MTase HVO_A0079 identified a number of homologs to the protein, including RM.Aco12261II (GenBank ADE57453), a Type IIG 6mA MTase in Aminobacterium colombiense DSM 12261 which has been identified via SMRT sequencing as targeting the motif CCRGA m6 G [32]. HVO_A0079 was also observed to share homology with RM.Fla104114II (GenBank BAV07385), a Type IIG 6mA MTase characterized via unpublished SMRT sequencing data on REBASE. Blastp also determined that HVO_A0079 shared homology with putative Type IIG 6mA MTases in Halorubrum californiensis DSM19288 (C463_0072; GenBank ELZ48543) and in Halophilic archaeon DL31 (RM.HarDL31ORF105P; GenBank AEN07377). These proteins range between 1117 to 1185 amino acids in length, which is similar to the 1088 amino acid length of HVO_A0079. These homologs are also all annotated on REBASE as Type IIG 6mA subtype α RM proteins, with the exception of C463_0072, which is not present in REBASE. A multiple sequence alignment of the amino acid sequences of these homologs (Figure 7) indicated significant sequence conservation in the central region of the alignment. Three large sections of this central region were observed via InterProScan to belong to the SCOP S-adenosyl-L-methionine-dependent methyltransferase superfamily domain SSF53335. InterProScan also identified three regions of the alignment which belong to the PRINTS adenine-specific DNA MTase signature domains PR00507, as well as a PFAM Eco57I domain PF07669, a domain observed in well-characterized Type IIG RM protein Eco57I [64]. A closer analysis of the alignment revealed the presence of FGG (conserved as AGG) and DPPY (conserved as NPPY) signature motifs in the order N-FGG-NPPY-C, which would follow the motif order N-FGG-TRD-DPPY-C observed in subtype α MTases according to REBASE [10]. However, no significant similarity was observed in the N-terminal region of these proteins, which is where the restriction endonuclease domain is typically located in Type IIG RM proteins. These results overall support the annotation of HVO_A0079 as a Type II 6mA subtype α MTase. InterProScan also identified three regions of the alignment which belong to the PRINTS adenine-specific DNA MTase signature domains PR00507, as well as a PFAM Eco57I domain PF07669, a domain observed in well-characterized Type IIG RM protein Eco57I [64]. A closer analysis of the alignment revealed the presence of FGG (conserved as AGG) and DPPY (conserved as NPPY) signature motifs in the order N-FGG-NPPY-C, which would follow the motif order N-FGG-TRD-DPPY-C observed in subtype α MTases according to REBASE [10]. However, no significant similarity was observed in the N-terminal region of these proteins, which is where the restriction endonuclease domain is typically located in Type IIG RM proteins. These results overall support the annotation of HVO_A0079 as a Type II 6mA subtype α MTase.

Deletion of HVO_0794, HVO_A0006, and HVO_A0237 Eliminates 4mC Methylation and Does Not Effect 6mA Methylation
In order to better understand the roles of HVO_0794, HVO_A0006, and HVO_A0237 in DNA methylation, the three genes were deleted in mrr deletion strain H1206, producing a triple deletion mutant (∆HVO_0794 ∆HVO_A0006 ∆HVO_A0237). The genome of this deletion mutant was sequenced via SMRT sequencing to determine the methylome and the results are listed in Table 5. In this strain, the C m4 TAG motif that is modified in the parental strain H26 [37] is no longer detected as methylated. Also, the 6mA motif GCA m6 BN 6 VTGC is modified in the triple deletion mutant, with 100% of the 410 motifs in the genome identified as methylated. Between studies, there was also a difference in the percent of motifs detected as methylated in H26 compared to the triple deletion mutant. In H26, only 316 of the 410 GCABN 6 VTGC motifs (~77%) were detected as methylated in [37], whereas in this study all 410 motifs are modified in ∆HVO_0794 ∆HVO_A0006 ∆HVO_A0237. This discrepancy is likely the result of a difference in sequence coverage, since the mean motif coverage and QV scores (confidence scores) were greater in the triple deletion mutant compared to H26. In ∆HVO_0794 ∆HVO_A0006 ∆HVO_A0237, the mean motif coverage for GCA m6 BGN 5 VTGC was 130.4, a~325% increase from the mean motif coverage of 30.7 in H26. The mean modification QV score for the 6mA motif was 213.0 in the triple deletion mutant, an increase of~274% from the H26 mean QV score of 57.0. Therefore, it is likely that the motifs were methylated completely in both strains, but that some of those motifs were not detected as modified in H26 due to the lower coverage and mean QV scores. Overall, these results indicate that deletion of HVO_0794, HVO_A0006, and HVO_A0237 abolishes methylation of C m4 TAG, but has no effect on methylation of GCA m6 BN 6 VTGC.

Deletion of the rmeRMS Operon Abolishes 6mA Methylation
The putative Type I operon rmeRMS was deleted in H. volcanii H1206, and sequenced via SMRT sequencing, in order to determine the role of the operon in DNA methylation. The results of the SMRT analysis for this strain (∆rmeRMS) are listed in Table 5. In ∆rmeRMS, the 6mA motif GCA m6 BN 6 VTGC is not detected as modified as it is in H26, and no other 6mA methylation is present [37]. Modification of C m4 TAG is still detected in the deletion strain. In ∆rmeRMS, 1199 of the 1342 CTAG motifs in the genome (~89%) are detected as methylated, These results are better than in Ouellette et al. [37] due to the better sequence coverage in ∆rmeRMS compared to H26, thus providing better detection of the methylated motifs. The mean motif coverage for C m4 TAG in this deletion mutant was 113.0. The mean modification QV score for the 4mC motif in ∆rmeRMS was 104.1. Overall, these results indicate that deletion of the rmeRMS operon eliminates methylation of the GCA m6 BN 6 VTGC motif.

Multi-RM Deletion Eliminates Detection of All DNA Methylation
A multi-RM deletion mutant, with all putative RM genes except for HVO_C0040 deleted from the strain, was also analyzed using SMRT sequencing to determine if the deletion of these genes resulted in elimination of DNA methylation in H. volcanii. The results of the SMRT analysis for this strain (∆RM) are listed in Table 5. The 6mA motif GCA m6 BN 6 VTGC identified in H26 is not detected as modified in this strain [37]. Also, the 4mC motif C m4 TAG is also not detected as methylated. No other motifs are detected as modified in this strain. These results indicate that all DNA methylation that can be detected by SMRT sequencing has been eliminated in the multi-RM deletion mutant. Although the remaining RM gene in this strain (HVO_C0040) encodes a MTase predicted to perform 5mC methylation which is difficult to detect via SMRT sequencing without Tet treatment [65], motifs of this type of methylation can still be weakly detected without Tet treatment. Since even weak detection of motifs was not observed in this strain, the results indicate that HVO_C0040 is not active as an MTase.

No Defect in Growth Occurs in the Multi-RM Deletion Compared to the Parental Strain
In E. coli dam − mutants, the lack of methylation results in growth defects compared to the wild-type strain [66]. In order to determine if the lack of RM genes resulted in a deficiency of growth in the ∆RM strain compared to the H26 parental strain, both strains were grown in Hv-YPC medium ( Figure 8). The results indicate that no significant difference in growth. Both strains entered log phase at 6 h, and although ∆RM initially had a slightly higher OD 620 when it entered log phase, this difference disappeared after~20 h of growth, and both cultures reached stationary phase at~36 h with similar OD 620 values (Supplementary Figure S1). The final OD 620 at for H26 after 72 h was 0.335, whereas for ∆RM the final OD 620 was 0.334. The difference between these two averages was not significant based on the standard error values and analysis of variance (ANOVA) single factor statistical analysis.
Overall, these results indicate that there is no detectable defect in growth in the ∆RM strain compared to the H26 strain.
Genes 2018, 9, x FOR PEER REVIEW 17 of 23 strain (∆RM) are listed in Table 5. The 6mA motif GCA m6 BN6VTGC identified in H26 is not detected as modified in this strain [37]. Also, the 4mC motif C m4 TAG is also not detected as methylated. No other motifs are detected as modified in this strain. These results indicate that all DNA methylation that can be detected by SMRT sequencing has been eliminated in the multi-RM deletion mutant.
Although the remaining RM gene in this strain (HVO_C0040) encodes a MTase predicted to perform 5mC methylation which is difficult to detect via SMRT sequencing without Tet treatment [65], motifs of this type of methylation can still be weakly detected without Tet treatment. Since even weak detection of motifs was not observed in this strain, the results indicate that HVO_C0040 is not active as an MTase.

No Defect in Growth Occurs in the Multi-RM Deletion Compared to the Parental Strain
In E. coli dam − mutants, the lack of methylation results in growth defects compared to the wildtype strain [66]. In order to determine if the lack of RM genes resulted in a deficiency of growth in the ΔRM strain compared to the H26 parental strain, both strains were grown in Hv-YPC medium ( Figure 8). The results indicate that no significant difference in growth. Both strains entered log phase at ~6 hours, and although ΔRM initially had a slightly higher OD620 when it entered log phase, this difference disappeared after ~20 hours of growth, and both cultures reached stationary phase at ~36 hours with similar OD620 values (Supplementary Figure S1). The final OD620 at for H26 after 72 hours was 0.335, whereas for ΔRM the final OD620 was 0.334. The difference between these two averages was not significant based on the standard error values and analysis of variance (ANOVA) single factor statistical analysis. Overall, these results indicate that there is no detectable defect in growth in the ΔRM strain compared to the H26 strain.

Discussion
In a previous study on DNA methylation in H. volcanii H26 [37], two motifs were identified as modified throughout the genome of the organism: the 4mC motif C m4 TAG and the 6mA motif

Discussion
In a previous study on DNA methylation in H. volcanii H26 [37], two motifs were identified as modified throughout the genome of the organism: the 4mC motif C m4 TAG and the 6mA motif GCA m6 BN 6 VTGC. These motifs were predicted to be methylated by a putative Type II 4mC MTase encoded by HVO_0794 and a putative 6mA MTase belonging to a Type I RM system encoded by the operon HVO_2269-2271 (rmeRMS), respectively. However, there are several annotated RM genes with no predicted motif recognition. In this follow-up study, we demonstrated through successive deletions of annotated RM genes that the C m4 TAG motif is methylated by the Type II MTase HVO_0794; the GCA m6 BN 6 VTGC motif is methylated by the Type I RM system RmeRMS; and that the other annotated MTases do not methylate under the conditions tested.
In mutants with HVO_0794 deleted from the genome, the SMRT sequencing analyses did not detect methylation of C m4 TAG or any other type of 4mC methylation, indicating that the MTase encoded by this gene is responsible for CTAG methylation since removal of this gene abolishes methylation of the motif. This result confirms predictions from previous studies [47,67] that suggested HVO_0794 is a CTAG MTase. Our bioinformatics analysis also supports the identification of this MTase as responsible for 4mC methylation, since the amino acid sequence has high similarity to previously characterized Type II 4mC CTAG MTases such as M.MthZI [33]. Although HVO_A0006 and HVO_A0237 were also deleted in the same strain as HVO_0794, neither of these genes have high similarity to 4mC CTAG MTases, and deletion of HVO_A0006 in a previous study [37] indicated that it does not affect cytosine methylation, ruling out these genes as candidates for C m4 TAG methylation. Based on our search of REBASE and NCBI, no cognate REase is encoded in the genome of H. volcanii, and deletion of the gene was not lethal as would be expected if there was a cognate REase, suggesting that HVO_0794 is an orphan MTase. The observation of an orphan CTAG MTase in H. volcanii was not unexpected based on previous work by Blow et al. [32], who found that predicted Type II CTAG orphan MTase gene families are common in the Halobacteria, occurring in 78% of halobacterial species.
Several of the halobacterial species examined by Blow et al. [32] which contained the CTAG MTase family also had a high CTAG motif density at their origins of replication, suggesting that this gene family may play a role in regulating DNA replication in the Halobacteria. Our analysis of H. volcanii showed a higher CTAG motif density surrounding oriC2, but not in regions near the other two origins; however, the oriC2 region also showed the enrichment of GATC motifs was even more pronounced (see Figure 3). It remains to be established, if the stretches with higher CTAG and GATC motif density have a selected function in H. volcanii, or if they reflect gene acquisition from a donor with different compositional bias.
H. volcanii has a lower GATC to CTAG ratio than organisms with homologs to the E.coli DNA adenine methyltransferase (Dam) which recognizes the GATC motif, and aids DNA repair via a methyl-directed mismatch repair system [21]. However, the decrease in the GATC to CTAG ratio in H. volcanii is due exclusively due to a drop in the frequency of the GATC motif, and not to an increase in the CTAG frequency. The role of CTAG methylation may not be of major importance for H. volcanii, as no growth defect was observed to occur in the ∆RM strain compared to the parental H26 strain. Considering that H. volcanii does not require origins of replication in order to grow efficiently [54], it is not too surprising that eliminating the putative role of HVO_0794 in regulating the origins of replication does not affect growth.
The absence of GCA m6 BN 6 VTGC methylation in deletion mutants without the rmeRMS operon indicated that these genes are responsible for 6mA methylation in H. volcanii. Our bioinformatics analysis also indicates that RmeRMS is a Type I RM system, since both the MTase subunit RmeM and specificity subunit RmeS are homologous to well-characterized Type I 6mA MTases and specificity subunits such as M.EcoKI and S.EcoKI [56]. The motif GCA m6 BN 6 VTGC resembles the type of sequences targeted by Type I systems, which are typically bipartite sequences with a gap of unspecified nucleotides in the middle [15]. Therefore, the observation that rmeRMS is a Type I RM system supports the identification of this operon as responsible for 6mA methylation in H. volcanii. Previous work by Ouellette et al. [37] had suggested that the RM gene HVO_A0006 might have a role in 6mA methylation, since SMRT sequencing of a HVO_A0006 deletion mutant identified an alteration in the 6mA motif (GCA m6 BGN 5 VTGC instead of GCA m6 BN 6 VTGC). However, our SMRT sequencing analysis of a deletion mutant without HVO_A0006 (∆HVO_0794 ∆HVO_A0006 ∆HVO_A0237) did not demonstrate any difference in the 6mA motif compared to the H26 parental strain. This difference is likely due to better sequence coverage in our data (~400x coverage for ∆HVO_0794 ∆HVO_A0006 ∆HVO_A0237 compared to~80x coverage for ∆HVO_A0006), allowing our analysis to identify more motifs as modified in the genome compared to the previous study [37]. This result, along with the observation that deletion of rmeRMS alone abolished detection of 6mA methylation, suggests that RmeRMS is solely responsible for adenine methylation in H. volcanii.
The presence of a restriction-subunit encoding gene (rmeR) indicates that the system can also cleave DNA at unmethylated target motifs, acting as a fully functional RM system. It is possible that this system functions in protecting H. volcanii from foreign DNA similar to RM systems in other organisms [68]. This defense system, in combination with clustered regularly interspaced short palindromic repeats (CRISPR-Cas system) [69,70], is likely advantageous to H. volcanii considering that haloarchaeoviruses are highly abundant in hypersaline environments [71,72]. RmeRMS may also be involved in regulating gene transfer, which in H. volcanii can occur within species as well as between species [73]. In E. coli, for example, the Type I RM system EcoKI has been demonstrated to reduce uptake via conjugation of unmethylated plasmids with EcoKI target sites [74]. A study by Lin et al. [75] indicated that RM systems could limit the size of DNA fragments that can recombine in Helicobacter pylori. Correlation between RM system occurrence and phylogenetic clusters was observed in Haemophilus influenzae, suggesting that RM systems are acting as barriers to genetic exchange between phylogenetic groups [76]. RM systems have also been hypothesized to drive population dynamics and diversification in Neisseria meningitidis [77]. Our results also indicate that homologs to RmeRMS, as well as the other predicted RM genes, do not occur as frequently in haloarchaeal species compared to the CTAG orphan MTase family genes; the RmeRMS system could possibly limit gene transfer that occurs with other individuals in the environment which lack the system, thus acting as a barrier to recombination for H. volcanii.
Our SMRT sequencing analyses indicate that deletion of HVO_0794 and rmeRMS is sufficient to eliminate detection of methylation in H. volcanii, indicating that the other predicted MTase genes in the organism (HVO_C0040, HVO_A0079, and HVO_A0237) do not contribute to methylation. The reason for the apparent inactivity of these genes is unclear, considering that our bioinformatics analyses indicate that these genes share homology with characterized MTase genes in other organisms. Inactive RM genes have been observed to occur in other organisms, such as those belonging to the MmeI RM gene family [78]. These inactive genes can be readily reactivated, and were hypothesized to exist in a population to confer a selective advantage to individuals when the population undergoes disruption from foreign parasitic DNA [78]. However, these MmeI RM genes were inactivated as a result of disruptive mutations which do not appear to be present in the predicted RM genes in H. volcanii. It is possible that these genes may still be active in H. volcanii but are only expressed under conditions not tested. However, a blastn search of the H. volcanii DS2 transcriptome data from Babski et al. [79] (sequence read archive (SRA) accession number SRP076059) using these three genes as queries suggested that they are expressed, although the search results do not indicate if functional protein products of these genes are produced. Interestingly, our results indicate that HVO_C0040 and HVO_A0079 are flanked by transposase genes similarly to HVO_A0237 [37]. Perhaps these genes are mobile genetic elements, as is the case with many RM genes [80], and they became non-functional when transferred into H. volcanii. Nevertheless, HVO_C0040, HVO_A0079, and HVO_A0237 do not appear to contribute to the methylome of H. volcanii under standard growing conditions. A possible exception to this list is HVO_C0040, the only remaining putative MTase gene in our ∆RM strain, which our bioinformatics analysis indicates is a 5mC MTase. Methylation patterns produced from 5mC MTases are typically difficult to detect with SMRT sequencing in the absence of Tet treatment [40,65].
However, 5mC methylation usually produces some signal via SMRT sequencing, yet we did not detect it in any of our multiple analyses including of the null mutant, leading us to think it is not methylating.
We also report in this study the construction of a MTase null mutant in H. volcanii. This strain (∆RM) has all putative RM genes deleted from the genome with the exception of HVO_C0040, and our SMRT sequencing analysis indicates that this strain has no genomic methylation. We anticipate that this strain will be useful for future studies that examine the impact of RM systems and DNA methylation on cellular processes in H. volcanii, in which the ∆RM strain can be compared to the parental strain H26 that has all the RM genes intact. This strain could also be useful for characterizing putative MTase genes in other halobacterial strains via gene knock-in and SMRT sequencing to determine the target sites for methylation. We expect that this strain will be a useful tool in the quest to better understand DNA methylation and RM systems in the Halobacteria and other archaeal organisms.
Supplementary Materials: The following are available online at www.mdpi.com/2073-4425/9/3/129/s1. Figure  S1. Growth curves of H26 and ∆RM when grown on Hv-YPC, represented by the average optical density (OD 620 ) readings of 24 cell culture replicates taken each hour for 72 h of growth.


extfieldsofstudy: Chemistry Medicine

s2fieldsofstudy: Biology Chemistry

text:
A single NaK channel conformation is not enough for non-selective ion conduction

NaK and other non-selective channels are able to conduct both sodium (Na+) and potassium (K+) with equally high efficiency. In contrast to previous crystallographic results, we show that the selectivity filter (SF) of NaK in native-like lipid membranes adopts two distinct conformations that are stabilized by either Na+ or K+ ions. The atomic differences of these conformations are resolved by solid-state NMR (ssNMR) spectroscopy and molecular dynamics (MD) simulations. Besides the canonical K+ permeation pathway, we identify a side entry ion-conduction pathway for Na+ permeation unique to NaK. Moreover, under otherwise identical conditions ssNMR spectra of the K+ selective NaK mutant (NaK2K) reveal only a single conformational state. Therefore, we propose that structural plasticity within the SF and the selection of these conformations by different ions are key molecular determinants for highly efficient conduction of different ions in non-selective cation channels.

S odium (Na + ) and potassium (K + ) are the most abundant monovalent cations in biological systems and play an essential role in the homeostasis and electrical activity of living cells. NaK belongs to the class of non-selective monovalent cation channels that conduct both K + and Na + ions with equal efficiency 1 . In contrast, K + selective channels are characterized by the remarkable ability to select K + over Na + by a ratio greater than 1000:1 2 . This high selectivity in K + channels is attributed to the narrow pore of the channel, known as the selectivity filter (SF). The conserved sequence TVGYG forms four K + -binding sites (S1-S4), which are composed from the backbone carbonyls and the threonine hydroxyl group 3 . The high conductance of K + selective channels has been explained in terms of either a multi-ion "knockon" mechanism 4-6 or the recently proposed "direct knock-on" mechanism 7 . Although NaK has a similar SF sequence (TVGDG) as K + selective channels, the mechanism by which NaK conducts both Na + and K + at extremely high rates remains under debate [8][9][10] .
High-resolution crystal structures of NaK revealed the conservation of only the S3 and S4 ion binding sites, with the upper part of the SF becoming a vestibule 11 . Notably, ion binding sites identical to those of NaK were observed in a recently solved cryo-electron microscopy structure of the biologically and physiologically important hyperpolarization-activated cyclic nucleotide-gated (HCN) channel 12 . HCN1 exhibits significant similarities to NaK in ion selectivity (K + /Na +~4 ), SF sequence, and structural architecture. These observations, together with mutagenesis studies on NaK, suggested that loss of contiguous binding sites (S1 and S2) is responsible for ion non-selectivity 13,14 . Moreover, no structural rearrangement was observed in the crystal structures of NaK with Na + , K + , or rubidium (Rb + ) and a single SF conformation was detected for NaK in bicelles at high temperature by solution-state NMR 15 . The different ions were therefore proposed to permeate through the same SF conformation with different preferential binding sites 11,16 .
In this study, we use a combination of ssNMR and advanced molecular dynamics (MD) simulations to investigate the detailed mechanism of ion non-selectivity in the NaK channel. Distinct from other methods in structural biology, ssNMR makes it possible to study membrane proteins in native-like lipid bilayers at room temperature and under physiological buffer conditions [17][18][19][20] . The MD simulations allow us to determine the conductive conformation of NaK for Na + and K + , and further enable us to investigate the detailed permeation mechanism of these ions [5][6][7] . Intriguingly, we identify two conformations of the SF in NaK, one of which is preferred by Na + and the other by K + . We further underline the functional importance of the conformational preferences induced by Na + or K + binding by using MD based permeation simulations. These results not only provide valuable insights into the conduction mechanism of NaK, but also have important implications for our understanding of the selectivity difference between the K + selective and non-selective channels.
Results ssNMR measurements and sequential assignments of NaK. We first prepared four uniformly [ 13 C, 15 N]-labeled NaK samples in liposomes: (A) at low ionic conditions (<1 μM) or with 50 mM of (B) Na + , (C) K + , or (D) Rb + . Single-channel conductance measurements verified that NaK in our lipid bilayer preparation is functional ( Supplementary Fig. 1). De novo sequential backbone and side chain chemical shift assignments were obtained from twodimensional (2D) and three-dimensional (3D) correlation spectra 21 of the Na + -containing sample ( Supplementary Figs. 2, 3). All visible signals in the spectra could be assigned to 56 residues, which are highlighted in the crystal structure and amino acid sequence of NaK (Fig. 1). Unassigned residues cluster to two different regions: (i) D66-S70 in the vestibule and extracellular entrance of the SF and (ii) the N-terminus up to F28 and from V91 until the C-terminus below the hinge region in the M1 and M2 helices. The signals from these regions cannot be detected in dipolar-coupling-based experiments, most probably due to conformational heterogeneity caused by strong structural dynamics. The signals from flexible residues, as observed in scalar-coupling based experiments, strongly correlate with the residue types found in the C-and N-termini ( Supplementary Fig. 4). However, the residues in the lower parts of the M1 and M2 helices are still not visible, implying that dynamics of these parts are relatively slow.
Two distinct states of the SF detected by ssNMR. The 2D 15 N i - 13 Cα i (NCA) and 15 N i -13 CO i-1 (NCO) correlation spectra of NaK at low ionic conditions showed two sets of signals for SF residues V64 and G65. Two sets of signals are also visible in the PDSD spectrum, highlighted for residue V64 in Supplementary Figure 5. The signal intensities suggest one dominant and one minor populated state, which are in very slow exchange compared to the μs-ms NMR timescale. For the upfield (lower chemical shift) nitrogen signals of V64 and G65, we observed ion-dependent chemical shift changes (low ionic vs K + vs Rb + ), while the downfield (higher chemical shift) nitrogen signals remained unaffected (Fig. 2). Since chemical shifts are highly sensitive to the local electronic environment induced by different ion types and occupancies in the SF 19,22 , this result strongly suggests that the upfield nitrogen signals correspond to a K + /Rb + -favored state and the downfield nitrogen signals correspond to an ion-free state.
The two detected states of the SF can be stabilized under different conditions. Low ionic conditions show the ion-free state as dominant, although the K + /Rb + -favored conformational state can still be observed. In contrast, this state becomes more dominant with increasing K + concentration (50 mM: Fig. 2, 150 mM: Supplementary Fig. 6). Intriguingly, ssNMR spectra of NaK in complex with Na + revealed a different equilibrium where the ion-free state is the major population even though ions are present. Furthermore, ssNMR spectra of the Na + sample revealed increased structural plasticity in the SF by significant line broadening for residue G65 (nitrogen dimension in the NCO spectrum) and by the appearance of additional peaks for V64 ( Fig. 2 and Supplementary Fig. 6). Under high Na + (150 mM) conditions, the signal intensity of the SF residues decreases further and both G65 and V64 show even more additional peaks that portray undefined conformations.
Structural characterization of the two SF states. The nitrogen chemical shift differences between the ion-free and K + /Rb + -favored states were up to 4 p.p.m. for V64 and G65. Such large nitrogen chemical shift differences are usually related to structural rearrangements associated with changes in hydrogen bonding. Carbonyl flipping in the SF should perturb the hydrogen bonding network and carbonyl flipped conformers have been observed in several crystal structures and simulations of KcsA 6,23,24 . In order to investigate alternative stable conformations of the SF, we undertook MD simulations on channels where we generated a set of backbone carbonyl flipped conformers for T62 to G65. Only the T62 flipped carbonyl conformer in all four subunits (Fig. 3c) was stable during the simulations under the ion-free condition (only ions to neutralize the system were present in the simulation box). This flipped conformation removes two hydrogen bonds for T63NH-T60CO and V64NH-V59CO, but permits the formation of a new hydrogen bond for V64NH-T62CO.
SsNMR distance measurements on NaK with K + provided strong support that the crystal structure and the flipped conformation correspond to the K + /Rb + -favored and ion-free states, respectively. Due to comparable rigidity in both conformations and similar neighboring amide nitrogen distances in the peptide chain, a 15 N-15 N PDSD spectrum with 5 s mixing time shows a clear cross-peak for T62′-T63′ with a similar intensity to the T62-T63 cross-peak ( Supplementary Fig. 7). In an NHHN spectrum with 200 μs proton-proton mixing time only close amide-proton contacts can be detected, because only distances below~4 Å give rise to strong cross-peaks 25 . For the K + /Rb + -favored state, we observed strong cross-peaks for all vicinal amide protons in the sequence 61 LTTVG 65 , which agrees perfectly with the crystal structure (Fig. 3). For the flipped conformation, only the distance between T62 and T63 is~4.5 Å and indeed lacks High-resolution ssNMR spectra of NaK. a SF in the crystal structure of KcsA (PDB ID: 1K4C) and NaK (PDB ID: 3E8H). b 2D NCA correlation spectrum of NaK in the presence of 50 mM K + with the signal assignment. SF residues are indicated in bold and the ion-free state is labeled with a prime (′) and in cyan to differentiate it from the K + /Rb + -favored state (no prime; in black). c Crystal structure of the open state of NaK displaying two opposing subunits. The hinge region is marked with an arrow and assigned residues are indicated on the structure in red. d Amino acid sequence of NaK with secondary structure elements based on the crystal structure indicated above the sequence. Assigned residues are shown in red and SF residues are accentuated by cyan arrows b c d 44   Single SF conformation observed in NaK2K. The different conformational states identified in NaK were not observed in the K + selective mutant NaK2K 13 (i.e., NaK D66Y and N68D double mutant) prepared under identical conditions. All SF residues could be assigned to a single signal, except T62 which shows an unconfirmed satellite peak ( Supplementary Fig. 8). NaK2K is characterized by a SF that is almost indistinguishable from other K + selective channels, such as KcsA, that contain four sequential ion binding sites. The upper part of the SF (Y66-S70), which cannot be detected for NaK, is clearly visible in the spectrum of NaK2K ( Supplementary Fig. 8). This indicates an increased rigidity for these residues, consistent with the formation of a canonical K + selective SF.
K + and Na + permeation by computational electrophysiology.
To check whether the different conformational states stabilized by K + or Na + are functionally relevant to the conduction properties of NaK, we performed atomistic MD simulations with a computational electrophysiology setup previously used to simulate ion permeation in K + selective channels (Fig. 4a) 7,26 . First, we conducted simulations using the crystal structure of NaK (PDB ID: 3E83) with KCl (Methods). Within 10 μs, we observed 18 inward and 19 outward permeation events at a transmembrane voltage of ±460 (±70) mV (Fig. 4b), indicating that the crystal structure is conductive for K + . Next, we employed the same setup to simulate Na + conduction in NaK, however no permeation could be detected within 10 μs. In accordance with previous simulations [27][28][29] , Na + ions remained tightly bound in the planes formed by backbone carbonyls of T63 or V64 for the duration of the simulations. To determine the conductive conformation of NaK for Na + , we performed a series of simulations using combinations of subunits in the flipped and crystal conformations (Supplementary Table 1). Intriguingly, only with mixed conformations were appreciable permeation rates for Na + detected that are comparable to K + conduction; for example, one flipped and three crystal subunits showed 23 inward and 9 outward permeations within 10 μs at a transmembrane voltage of ±560 (±40) mV (Fig. 4b). Although the simulated conductance (0.6 pS for inward K + , 0.8 pS for inward Na + ) is about 44 times lower than observed experimentally (35 pS) 13 , the conductance ratio between the previously simulated MthK channel (MthK: 7.2 pS for outward K + ) 7 and NaK from our study is about 9. This agrees relatively well with the experimentally observed singlechannel conductance ratio (MthK/NaK: 6.3) 13,30 .
Asymmetric SF revealed during Na + conduction. The SF shows a break in fourfold symmetry by having one or multiple subunits in the T62CO flipped conformation when conducting Na + . This leads to an asymmetric deformation of the SF during the simulations, as reflected by the distribution of the distance between opposing backbone carbonyls for residues T63 and V64 (Fig. 5). Ion occupancy ( Supplementary Fig. 9) as well as 2D potential of mean force ( Supplementary Fig. 10) were calculated from the permeation simulations, which had either a mixture of flipped and crystal conformations or the crystal conformation in all subunits. The comparison suggests that the structural deformation in the SF caused by T62CO flipping significantly reduces the energy barriers for Na + conduction at sites B 23 and B 34 .
Side entry ion-conduction pathway for Na + permeation. Surprisingly, we found two distinct ion-conduction pathways for Na + and K + permeation through the SF in the MD simulations. Outward and inward K + currents (see Supplementary Fig. 11 for definition) are dominated by the canonical K + conduction pathway (Supplementary Movies 1 and 2), where multiple K + ions pass through the entire narrow pore of the SF in single file. However, for Na + conduction in both directions, ions use side entries instead of passing through the vestibule and pore entrance ( Fig. 4b and Supplementary Movies 3 and 4). Notably, the additional binding site at the side entry (S side ) (Fig. 4c) is composed of two backbone carbonyls from D66 and F69 and the side chain oxygen from either S70 of the same or D66 of the neighboring subunit. Therefore the side entry pathway is a result of the perpendicular reorientation of the backbone carbonyl from D66 in NaK, which is a tyrosine residue in K + selective channels that points toward the fourfold symmetry axis of the channel. Ion occupancy and ion hydration states. During the outward and inward conduction of Na + and K + , there were significant differences in the ion occupancies as well as hydration states of ions in the SF. During Na + conduction, only one fully hydrated ion occupied the SF most of the time, while one or more Na + ions bound at S side . The ion at B 23 is weakly populated and short lived ( Supplementary Fig. 9). Water molecules were co-transported with the cations during these simulations. In contrast, during K + conduction, usually two dehydrated K + ions simultaneously occupied the S3 and S4 binding sites together with one fully hydrated ion in the vestibule (Supplementary Fig. 12). Additionally, we noticed that a few water molecules were able to be cotransported during the simulations of K + conduction in NaK ( Supplementary Movies 1 and 2).

Discussion
NaK is a prokaryotic non-selective cation channel that is able to conduct both Na + and K + with equally high efficiency. Due to the significant similarity in their SF sequences as well as their non-selective nature, NaK has been considered a bacterial homolog of eukaryotic cyclic-nucleotide-gated channels. In the present study, a combined use of ssNMR and MD simulations revealed two different conformational states of the SF that exist simultaneously in NaK and are stabilized by different ions. The K + /Rb + -favored state was confirmed as the crystal conformation, while the ion-free state corresponds to the T62CO flipped conformation. Most intriguingly, Na + favors the ion-free state and results in an asymmetric deformation of the SF.
The structural complexity of Na + binding was experimentally supported by the appearance of additional peaks in ssNMR spectra of NaK at high Na + concentration (150 mM NaCl, Supplementary Fig. 6), which were not visible with K + at the same concentration. This indicates an importance of conformational dynamics during Na + conduction, since the additional peaks appear only under high Na + conditions. We compared the distance information of the ion-free and K + /Rb + -favored states (obtained from the ssNMR spectra) with the SF conformations seen in KcsA. The conductive and collapsed conformations of KcsA show no significant difference in the distances between neighboring amide protons ( Supplementary  Fig. 13). Thus, the conformational dynamics observed in the SF of NaK are absent in the K + selective channel KcsA. However, comparison of the dipolar-coupling and scalar-coupling based NMR experiments suggests that the dynamics of the C-and Nterminal regions are relatively slow in NaK. This behavior is most likely important for channel gating, which is in accordance with  Fig. 4 MD simulations of K + and Na + permeation. a The simulation system consists of two membranes with NaK (PDB ID: 3E83), surrounded by water and ions, and exhibits a transmembrane voltage gradient (see also Supplementary Fig. 11). b Inward and outward permeation events of K + (red) and Na + (blue) vs time. K + simulations were performed with the crystal conformation, while Na + simulations were conducted using a mixture of crystal and T62CO flipped conformations in different subunits. Ion permeation profiles for outward K + and inward Na + are shown on the right. c S side coordination of Na + during the inward Na + simulations. the activation gating mechanism of KcsA by allosteric transmembrane coupling 31,32 .
To further validate the functional relevance of the conformational plasticity of the SF in NaK, we performed the ssNMR experiments also on a K + selective NaK double mutant (NaK2K) under otherwise identical conditions. Here we observed only a single conformational state of the SF, which is in strong contrast to NaK. Furthermore, the upper part of the SF (Y66-S70) became visible, as would be expected from a KcsA-like SF conformation. These results eliminate the possibility that the observed conformational heterogeneity in NaK is a result of our sample preparation.
Computational electrophysiology enabled the determination of the conductive conformation for Na + , while the crystal structure is only conductive for K + . Efficient Na + conduction in NaK requires conformational plasticity and the resulting asymmetric deformation of the SF.
Most importantly, we show that the permeation mechanisms of K + and Na + are remarkably different. In NaK several K + ions move in a single file along the axis of the SF via a loosely coupled direct knock-on mechanism. Although a few water molecules could be co-transported during K + conduction, the permeation mechanism is similar to that proposed for K + selective channels using the same computational approach 7 . Intriguingly, for the inward and outward Na + current, we identified a different pathway in which Na + uses a peculiar side entry into the SF (Fig. 6c). The additional binding site at the side entry is attributed to the perpendicular reorientation of the backbone carbonyl of D66 in NaK. We also observed the side entry pathway for K + permeations in a set of simulations, where V64CO flipping resulted in a similar kind of asymmetric deformation of the SF. Nevertheless, since the V64CO flipped conformation in the presence of K + was not supported by the ssNMR data, we do not consider this side entry pathway to be physiologically relevant for K + .
Moreover, in marked contrast to the K + conductive state, the lower part of the SF in the Na + conductive state becomes wider (Fig. 5). The wider pore enables only one fully hydrated Na + ion to reside in the SF with a short occupation time at B 23 . This result further indicates that during Na + permeation ions do not need to be fully dehydrated. The Na + permeation mechanism in NaK resembles the one that was observed for Na + conduction through the open conformation of a bacterial voltage-gated Na + channel 33 . As the SF of the NavM channel 34,35 is significantly wider than that of K + selective channels, the ions are hydrated at all sites during the Na + permeation simulations. Previous studies suggested that dehydration of the smaller Na + ions requires much more energy than of larger K + ions. Conduction through the SF without full dehydration is thus beneficial for Na + ions 36,37 . Furthermore, it should be noted that the SF of the less selective Na + channels, as well as of several other recently structurally elucidated non-selective ion channels 38-40 , is considerably wider than the SF of highly K + selective channels. Consequently, we propose that the structural plasticity and the ability of adapting the size of the SF in NaK are key molecular determinants for its ion non-selectivity. Previous NMR and crystallographic studies on both NaK and KcsA did not reveal the same degree of conformational plasticity in the SF 15,17 . We believe our model complements the previously proposed thermodynamic and kinetic binding mechanisms 10 as well as the single-barrier mechanism of non-selectivity in NaK 11,13 and thereby provides a more thorough understanding of ion selectivity in K + selective and non-selective channels.

Methods
Sample preparation. The NaK gene without the first 17 amino acids was cloned into a pET28a expression vector (Supplementary Table 2) with an N-terminal hexahistidine (His 6 ) sequence. To prepare the [ 13 C, 15 N]-labeled samples, the protein was expressed in BL21(DE3)-Gold cells (Agilent) in M9 medium with 13 Cglucose as the sole carbon source and 15 NH 4 Cl as the sole nitrogen source. The cultures were grown for 12 h at 20°C after induction. The membrane fraction was isolated and solubilized in 20 mM DM (n-decyl-beta-D-maltopyranoside). The protein samples were purified by Ni-NTA affinity chromatography and size exclusion chromatography either without cations or in the presence of 50 or 150 mM Na + , 50 or 150 mM K + , or 50 mM Rb + . The concentration of purified NaK was determined by the A280 method with the calculated extinction coefficient. The NaK tetramers were reconstituted into asolection liposomes at a lipid/tetramer ratio of 40:1 using dialysis. Proteoliposomes were collected by ultracentrifugation at 300,000×g for two hours at 4°C.
Single-channel conductance measurements. After dialysis, the protein/lipid mixture was extruded 20 times through a 100 nm polycarbonate membrane filter (Whatman, Newton, MA) with an Avanti mini-extruder to make unilamellar liposome vesicles. Channel conductance measurements in planar lipid bilayers were conducted on an Ionovation Compact device (Osnabruck, Germany). Two polycarbonate compartments with a volume of 1.2 mL were separated by a sheet of TEFLON foil with 25 μm thickness and 50−100 μm aperture diameter. The artificial lipid bilayer was formed in the aperture through the method of painting. The cis-chamber was filled with 15 mM NaCl, 150 mM KCl, 5 mM Mops, 5 mM Tris-HCl, pH 7.0, and the trans-chamber with 150 mM NaCl, 15 mM KCl, 5 mM Mops, 5 mM Tris-HCl, pH 7.0. Voltage was applied across the bilayer using Ag/AgCl electrodes immersed in each chamber. With the grounded cis-compartment, a positive potential indicating a higher potential in the trans-chamber was applied. Planar lipid bilayer formation was monitored optically or by capacitance measurements. After successful formation of a stable bilayer in the aperture, NaK proteoliposomes were added to the trans-chamber next to the bilayer. Fusion of the protein and the planar lipid bilayer was detected through observation of channel conductance. Signal acquisition and analysis were performed using the pCLAMP software (Axon Instruments). Single-channel conductance events were identified automatically and analyzed using the Clampfit 10.0 software (Axon Instruments  Fig. 6 Schematic model of the non-selective ion-conduction mechanism in NaK. a Equilibrium of NaK conformations without ions or in the presence of Na + or K + . Plasticity in the SF enables Na + binding, further resulting in an asymmetric deformation of the SF. b K + conduction pathway from the MD simulations. c Na + conduction in both directions uses the four side entries (dashed line). The additional binding sites (S side ) at the side entry are labeled as S S zero, indicating that NaK could conduct both Na + and K + ions with low ion selectivity ( Supplementary Fig. 1).
ssNMR spectroscopy and analysis. The proteoliposome pellets were transferred into 3.2 mm magic-angle spinning (MAS) rotors and used for ssNMR analysis. SsNMR spectra required for chemical shift assignment were recorded on a 16.4 T wide-bore NMR spectrometer (700 MHz, Bruker BioSpin) equipped with a 3.2 mm triple-resonance Efree MAS probe. Spectra were recorded at 17 kHz MAS rate and calibrated with external DSS (4,4-dimethyl-4-silapentane-1-sulfonic sodium salt) as a proton chemical shift reference. Intra-residue spin system correlations were obtained from 3D NCACB and band-selective homonuclear cross-polarization (BSH-CP) 41,42 based NCACO spectra, whereas inter-residue connections were provided by 3D CANCO and BSH-CP based NCOCA and NCOCACB spectra. As an illustration, a sequential walk from A53 to S57 is shown in Supplementary Fig. 2. The inter-residue Cα i -Cα i-1 connections were further confirmed by a 2D CANCOCA spectrum. Subsequently, the side chain carbon chemical shifts were assigned with a 3D NCACX spectrum. A 13 C-13 C proton-driven spin diffusion (PDSD) 43 spectrum was recorded at 11 kHz MAS with 20 ms mixing time to assist the side chain assignment (Supplementary Fig. 3). An INEPT-based 1 H-13 C correlation spectrum was recorded at 11 kHz MAS to detect flexible residues. The complete chemical shift list for NaK in the presence of 50 mM Na + is provided in Supplementary Table 3.
For the detection of short (<4 Å) distances between amide protons in the SF of NaK, an NHHN spectrum was recorded at 16.4 T and 11 kHz MAS using 200 μs proton-proton mixing time. As a reference, a 2D 15 N-15 N PDSD spectrum was recorded at 11 kHz MAS with 5 s mixing time.
The NCA and NCO spectra that were used to detect cation-specific chemical shift changes in NaK were recorded on an 18.8 T wide-bore NMR spectrometer (800 MHz, Bruker BioSpin) equipped with a 3.2 mm triple-resonance Efree MAS probe at 11 kHz MAS rate. The resulting spectra are shown in main text Fig. 2 (excerpts) and Supplementary Fig. 6 (full spectra). The chemical shifts of the SF residues under different ionic conditions (<1 μM; 50 mM Na + ; 50 mM K + ; 50 mM Rb + ) are listed in Supplementary Table 4.
During all experiments, the sample temperature was kept constant at 4-8°C, as measured by the temperature-dependent position of the water resonance 44 . Highpower proton decoupling using the sequence SPINAL-64 45 with a radio frequency amplitude of about 83 kHz was applied during evolution and detection periods.

Molecular dynamics simulations.
In all simulations, the NaK channel was embedded in a hydrated 1-palmitoyl 2-oleoyl-phosphatidylcholine (POPC) lipid bilayer with all atoms, including those in the lipids and water molecules, represented explicitly 46 .
The starting structure in the simulations was based on the high-resolution crystal structure of NaK (PDB ID: 3E83). All titratable residues were left in the dominant protonation state at pH 7.0. For the computational electrophysiology simulations, the F92A mutant instead of the wild-type (WT) channel was employed, as previously performed electrophysiology experiments revealed a significant increase in ion-conduction rates for the F92A mutant as compared to the WT channel 47 .
MD simulations were carried out using GROMACS 5.0 48,49 . We used the AMBER99sb force field 50 and the SPC/E water model 51 for the equilibrium and production simulations of the proteins. Parameters for ions and lipids were derived from references 52,53 . Short-range electrostatic interactions were calculated with a cutoff of 1.0 nm, whereas long-range electrostatic interactions were treated by the particle mesh Ewald method 54,55 . The cutoff for van-der-Waals interactions was set to 1.0 nm. Virtual sites were used for protein hydrogen atoms during the simulations, which allows for a 4 fs time step for integration. The simulations were performed at 300 K with a velocity rescaling thermostat 56 . The pressure was kept constant at 1 bar by means of a semi-isotropic Parrinello-Rahman barostat 57,58 . All bonds were constrained with the LINCS algorithm 59 .
To equilibrate the system, a short 20 ns simulation was performed by positionrestraining all heavy atoms of the NaK channel with a force constant of 1000 kJ mol −1 nm −2 to the starting structure. The system was further equilibrated for 20 ns by releasing the position restraints.
SsNMR spectra revealed two well-separated signals for SF residues V64 and G65, suggesting that the two corresponding conformations of the SF are in very slow exchange on the NMR timescale. This result further indicates that interconversion from one conformation to the other is not expected to take place within ns-μs MD simulations. As flipping of backbone carbonyls has been observed in some structures of mutants and in simulations of KcsA 6,23,24 , we assumed that the second signal set of the ssNMR spectra corresponds to a backbone carbonyl flipped conformer. In order to investigate which carbonyl flipped conformers are thermodynamically stable, we first generated a set of carbonyl flipped conformers by manually flipping the backbone carbonyls of T62, T63, V64, and G65 in the SF by~180°using the program Coot 60 . Two other residues in the SF (D66 and G67) were not tested, as no NMR signals were observed for these residues, indicating a fast conformational exchange on the ns-μs timescale. The manually flipped conformers were equilibrated first by positionrestraining all heavy atoms of the SF with a force constant of 1000 kJ mol −1 nm −2 to the starting structure. Starting from these equilibrated structures, we further tested the stability of the carbonyl flipped conformers by performing 5 independent runs of 20 ns MD simulations without position restraints. Among all carbonyl flipped conformers, only the T62CO flipped conformer was stable during these simulations, while for the other conformers the backbone carbonyl quickly returned to the crystal conformation.
For the computational electrophysiology simulations 61 , the equilibrated system was duplicated along the z direction and transmembrane potential gradients were generated by introducing two more K + or Na + ions in compartment (a) than in compartment (b) (Supplementary Fig. 11). During the MD simulations, the number of ions in these two compartments was kept constant by an additional algorithm 61 . The charge difference leads to a positive voltage difference for the upper channel and a negative voltage difference for the lower channel ( Supplementary Fig. 11). The resulting membrane potential can be calculated by double integration of the charge distribution using the Poisson equation as implemented in the GROMACS tool g_potential 62 . During the simulations an outward permeation event was counted when an ion moved from the cavity to the SF and another ion left the SF, while an inward permeation event was counted when an ion moved from the extracellular side to the SF and another ion left the SF. Snapshots of the NaK simulations including ions and water molecules are shown in Supplementary Figure 12.
Data availability. Data supporting the findings of this manuscript are available from the corresponding authors upon reasonable request. The NaK chemical shifts have been deposited in the Biological Magnetic Resonance Data Bank under accession codes 27219 (the ion-free conformation) and 27220 (the K + /Rb + -favored conformation).