Genomic DNA was prepared using the Sigma GenElute bacterial genomic DNA kit according to the manufacturer's instructions. Oligonucleotide primers were purchased from Integrated DNA Technologies (Singapore). PCR was carried out using GoTaq polymerase (Promega, Madison, USA) or KOD HotStart DNA polymerase (Merck, Darmstadt, Germany) according to the manufacturer's instructions. Plasmid DNA was prepared using the Promega Wizard Mini Prep kit (Madison, USA). Restriction enzymes were purchased from New England Biolabs (NEB, Ipswich, MA, USA). DNA ligase, polynucleotide kinase (PNK) and Antarctic phosphatase used during cloning were all purchased from NEB (Ipswich, MA, USA). Purification of PCR products was carried out using the Qiagen PCR purification kit according to the manufacturer's instructions.

A. pleuropneumoniae isolates were grown on BBL™ blood agar base supplemented with thiamine HCl (0.0005%), NADH (0.0025%), oleic acid bovine albumin complex (5%) consisting of 4.75% bovine serum albumin (Fraction V) in normal saline (with the saline containing 5% NaOH and 0.06% oleic acid) and heat inactivated horse serum (1%) (termed BA/SN) (41) or brain heart infusion (BHI, Oxoid) agar containing β-nicotinamide adenine dinucleotide, NAD (0.01%) at 37°C for 24 h (42). Escherichia coli strains (DH5α and BL21) were grown overnight in Lysogeny broth (LB) supplemented with ampicillin at 100 μg/ml or chloramphenicol at 20 μg/ml as required at 37°C in a shaker at 200 rpm (39). For mutant selection, A. pleuropneumoniae strains were grown overnight on BHI-NAD plates with 1 μg/ml of chloramphenicol. The Type III methyltransferase genes were cloned into the pET15b vector (Novagen). The prototype strains of A. pleuropneumoniae used in this study are listed in Table 1.

A total of 210 previously sequenced genomes, and all publicly available genomes from the RefSeq database in NCBI GenBank were used to determine the distribution of Type I and Type III R-M systems in A. pleuropneumoniae (Supplementary Table S1). We used as search seeds the conserved hsdR and hsdM region of the Type I RM system from AP76 (GenBank locus tags APP7_1526-APP7_1531) (25), and prototype Type III modP1-4 and modQ genes from AP76 (GenBank locus tags APP7_0747), JL03 (GenBank locus tags APJL_0704) (Table 1). A tblastn (NCBI GenBank) search against a custom blast database of our collection of A. pleuropneumoniae was used to score the presence or absence of each system, where presence was determined by a >90% nucleotide identity over 90% of the length of the gene. To create a core genome alignment and phylogeny of our collection, we first annotated open reading frames using Prokka (v1.12) (43). Then we used the pangenome software Panaroo (v1.2.7) (44) with default parameters to create a core genome alignment. A neighbour joining tree was created from this alignment using the functions dist.dna and nj in the R package ape (45,46) The distribution of the Type I and III systems were annotated on to the phylogeny using the interactive Tree of Life (iTOL) (47).

An allele specific multiplex PCR was developed based on the TRD region of the prototype strains to identify different Type III methyltransferase genes and alleles present in un-sequenced isolates. A total of 250 field isolates of A. pleuropneumoniae submitted from 2000 to 2019 to the reference identification and serotyping service provided by the Microbiology Research Group (UQ, Brisbane, Australia) and serovar 1 to 15 A. pleuropneumoniae reference serovar strains were screened for Type III mod genes and their allelic variants (Supplementary Table S2) (48,49). Multiplex PCR was carried out using GoTaq DNA polymerase (Promega) according to manufacturer's instructions using the primers Mod-TRD-P1-F and Mod-TRD-P1-R, Mod-TRD-P2-F and Mod-TRD-P2-R, Mod-TRD-P3-F and Mod-TRD-P3-R, Mod-TRD-P4-F and Mod-TRD-P4-R, Mod-TRD-Q-F and Mod-TRD-Q-R, listed in Supplementary Table S3. Cycle conditions were as follows: initial denaturation at 95°C for 2 min, followed by 30 cycles of denaturation at 95°C for 30 s, annealing at 63°C for 30 s, extension at 72°C for 1 min, with a final extension at 72°C for 10 min. Samples were checked on 1.5% (w/v) agarose gels buffered with 1x Tris-borate-EDTA.

To quantify the ON-OFF ratio of each Type III mod gene, a mod specific primer pair was designed to amplify over the simple sequence repeat (SSR tract), using methodology as previously described (39). Briefly, a gene specific primer pair was designed for each of the two mod genes present (App-ModP-F-FAM and App-ModP-R or App-ModQ-F-Yakyel and App-ModQ-R, Supplementary Table S3) with the forward primer labelled with a fluorescent label to allow quantification using fragment length analysis (Life Technologies Australia Pty Ltd). Enrichment of strains expressing each Type III mod variant to generate populations of bacteria containing ∼90% ON or ∼90% OFF was carried out by single colony screening and enrichment (28). PCR was carried out using Go-Taq DNA polymerase (Promega) according to manufacturer's instructions with following cycling conditions: initial denaturation at 95°C for 2 min, followed by 30 cycles of denaturation at 95°C for 30 s, annealing at 52°C for 30 s, and extension at 72°C for 30 s, with a final extension at 72°C for 5 min. PCR fragments were sent to Australian Genome Research Facility (AGRF, Brisbane, Australia) and analysed with Peak Scanner Software v1.0 (Life Technologies Australia Pty Ltd).

Triplicate cultures of A. pleuropneumoniae enriched ModP and ModQ ON/OFF variants were grown to mid-log (OD₆₀₀ ∼ 0.4–0.7) in BHI broth supplemented with NAD (0.01%). RNA was prepared from 2 ml of bacterial culture using Trizol reagent (Thermo-Fisher) according to the manufacturer's instructions. This RNA was used to synthesize cDNA using Protoscript II Reverse Transcriptase (NEB) and random hexamers (NEB) according to the manufacturer's instructions. A reverse transcriptase reaction without Protoscript II was used as a negative control.

Semi-quantitative RT-PCR was carried out using GoTaq DNA polymerase (Promega) with 16S rRNA primers (16S F and 16S R) and one of four Type I allele- specific primer pairs (alleles A, B, C, D, E, F and G) (Supplementary Table S3) in a multiplex reaction using following cycling conditions: 95°C for 2 min, followed by 30 cycles at 95°C for 30 s, annealing at 54°C for 20 s, extension at 72°C for 45 s and final extension at 72°C for 5 min (39). Samples were checked on 1.5% (w/v) agarose gels buffered with 1× Tris–borate–EDTA.

Allelic variants of each hsdS gene in A. pleuropneumoniae strain AP76 and JL03 were quantified using FAM labelled PCR coupled to a restriction digestion, based on methods as previously described (34,39). Specifically, PCR amplification of the whole hsdS locus (3.7 kb) was carried out from extracted genomic DNA using GoTaq DNA polymerase (Promega) with primers App_T1_For_FAM and App_T1_Rev (Supplementary Table S3) according to manufacturer's instructions. PCR cycle conditions were as follows: initial denaturation at 95°C for 2 min, followed by 30 cycles of denaturation at 95°C for 30 s, annealing at 52°C for 30 s, and extension at 72°C for 4 min, with a final extension at 72°C for 5 min. The PCR product was then subjected to restriction digestion with BccI and NheI, with the digestions predicted to produce different sized fragments for each of the variant forms (Supplementary Figure S1). DNA fragments were sent to the Australian Genome Research Facility (AGRF, Brisbane, Australia) or the Griffith University DNA Sequencing Facility (GUDSF) and analysed with Peak Scanner Software v1.0 (Life Technologies Australia Pty Ltd). The area of the peak given by each labelled fragment, each corresponding to the prevalence of one of the variant forms, was quantified using Peak Scanner v1.0.

The modP (modP1 and modP2) and modQ genes were PCR amplified using primers ModP-oE-F and ModP-oE-R, ModQ-oE-F and ModQ-oE-R (Supplementary Table S3) with KOD Hot-start DNA polymerase (Merck) according to manufacturer's instructions. The resulting product was cloned into the NdeI and BamHI site of the pET15b vector in such a way that the protein was expressed with an N-terminal His-tag, and designated as pET15b::modP1, pET15b::modP2, and pET15b::modQ. After confirmation of the correct construct by colony PCR (primers: T7 forward and App-ModP-long for ModP and T7 forward and App-ModQ-long for Mod Q; Supplementary Table S3) and restriction digestion, proteins were over-expressed in E. coli BL21 at 37°C with Isopropyl β-D-1-thiogalactopyranoside (IPTG) induction (0.5 mM) with 200 rpm shaking. A no induction control was performed under identical conditions. The presence of over-expressed protein was confirmed using monoclonal anti-poly histidine alkaline phosphatase antibody (Sigma-Aldrich) at 1:20 000 dilution.

All Mod knockout stains were generated by replacing the TRD region with a chloramphenicol acetyltransferase gene (cat) (50), previously used to make insertional knockouts in A. pleuropneumoniae (51). The TRD region was first removed from the pET15b::modP1 and pET15b::modQ by inverse PCR using the primers ModP_TRD_remove_F and ModP_TRD_remove_R; ModQ_TRD_remove_F and ModQ_TRD_remove_R, (Supplementary Table S3) and the amplified cat gene (Cm_USS_For and Cm_USS_Rev, Supplementary Table S3) was inserted by blunt end ligation. The constructs created were linearized and transformed into representative A. pleuropneumoniae strains via MIV transformation (52). Successful transformants were confirmed by PCR (App-ModP-F-FAM and CM-F-check; App-ModQ-F-Yakyel and CM-F-check, Supplementary Table S3) and sequencing.

TRD-less Mod proteins were constructed via inverse PCR with primers ModP_TRD_remove_F and ModP_TRD_remove_R; ModQ_TRD_remove_F and ModQ_TRD_remove_R (Supplementary Table S3) using the vector pET15b::modP1 and pET15b::modQ as a template. The PCR reaction was carried out using KOD Hot Start DNA polymerase (Merck) according to manufacturer's instruction. The conserved 5′ and 3′ region was then fused and the TRD-less proteins over-expressed using E. coli BL21. Over-expression was carried out overnight at 37°C (200 rpm) using IPTG (0.5 mM). After over-expression, proteins were purified using Talon resin (Takara Bio) by resuspending the cell pellet in lysis buffer (50 mM phosphate buffer and 300 mM NaCl) containing 0.1 mg/ml lysozyme and 0.5% SDS and multiple rounds of sonication. Proteins were eluted from the resin using gradient imidazole concentration (10–250 mM) in 1× binding buffer. After analysing fractions by SDS-PAGE and Western blot, pure fractions were concentrated and buffer exchanged using centrifugal concentrators (Amicon). TRD-less Mod proteins for both ModP and ModQ were then further purified by electroelution (100 mA for 30 min) into 1× MOPS buffer using a mini whole gel eluter (Biorad, Gladesville, NSW, Australia). Protein fractions were collected from each well and analysed by silver staining (Sigma) and Western blot using monoclonal anti poly histidine alkaline phosphatase antibody (1:20 000 dilution). Pure Mod TRD-less proteins (ModP and ModQ) were dialysed overnight at 4°C in 1× PBS twice.

Anti ModP and ModQ sera (primary mouse antibody) were generated by immunizing five female BALB/c mice (6–8-week-old). A 50 μg aliquot of electroeluted purified proteins (TRD-less ModP and ModQ proteins) in alum were administered sub-cutaneously per mouse on day 0, 14, 28, 42. Terminal bleeds were collected on day 58. Serum was separated by centrifugation and stored at –20°C in 50% glycerol. These antibodies (anti-ModP and anti-ModQ sera) were used as primary antibody. All animal work was approved by Griffith University Animal Ethics Committee (GLY/16/19/AEC).

Whole cell lysates of each of the enriched Mod variants were prepared from mid-log (OD₆₀₀ ∼ 0.4-0.7) cultures grown in BHI-NAD broth. Samples were normalized to OD₆₀₀ of 2.0 and prepared for SDS-PAGE using Novex Bolt™ LDS sample buffer and β-mercaptoethanol (final concentration 10% [v/v]). Samples were boiled for 20 min at 95°C and centrifuged at 14 000 × g for 20 min. Samples were run at 150 V for 1 hr on Novex Bolt™ Bis–Tris Plus precast Gels (4–12%) with 1× MOPS buffer (Life Technologies Australia Pty Ltd) and transferred to a nitrocellulose membrane (Bio-Rad, California, United States) at 15 V for 90 min. The membrane was then blocked overnight with 10% (w/v) skim milk in 1× Tris-buffered saline with Tween 20 (0.05%) (TBST) with shaking. The membrane was washed with 1× TBST 3–5 times for a total of 20 min, then incubated with primary mouse antibody raised against TRD-less ModP and ModQ proteins at 1:1000 dilution for 1 h at room temperature. The membrane was washed 3–5 times with 1× TBST for a total of 20 min and then incubated for an hour with 1:3000 dilution of secondary antibody (goat anti-mouse alkaline phosphatase conjugate, Sigma). The membrane was washed again with 1× TBST 3–5 times for a total of 20 min and developed using 5-bromo-4-chloro-3-indolylphosphate (BCIP)/nitro-blue tetrazolium (NBT) (Roche) according to the manufacturer's instructions.

SMRT sequencing and methylome analysis was performed as described previously (53,54) using genomic DNA isolated using the Sigma GenElute kit, or plasmid DNA isolated with the Qiagen miniprep kit. DNA was sheared (plasmid: ∼0.5–2 kb; gDNA: 5–10 kb) using g-TUBEs (Covaris, Woburn, MA, USA) and sheared DNA was used to prepare SMRTbell template sequencing libraries. After end repairing, DNA was ligated to hairpin adapters and a combination of Exonuclease III (NEB; Ipswich, MA, USA) and Exonuclease VII (USB; Cleveland, USA) were used to degrade incompletely formed SMRTbell templates. Primer was annealed and samples were sequenced for long insert libraries on the PacBio Sequel system. SMRT sequencing and methylome analysis was carried out at SNPSaurus (University of Oregon, USA).

Genomic and plasmid DNA was isolated as above for SMRT sequencing. Nanopore sequencing libraries were prepared from 400 ng of DNA per sample with the Rapid DNA Barcoding kit (RBK004) according to manufacturer's instructions. A 500 ng pool of 6 indexed plasmid libraries was sequenced on 1 PromethION R9.4.1 flow cells (FLO-PRO002) on a PromethION 24 instrument running MinKNOW 22.03.4. Reads were base called in super-accuracy mode (dna_r9.4.1_450bps_sup_prom) with demultiplexing, using Guppy 6.0.7 and default parameters. Read depth was downsampled to 400× per plasmid with ont_fast5_api v4.0.1 (https://github.com/nanoporetech/ont_fast5_api), first filtering for read quality with Filtong v0.2.1 (https://github.com/rrwick/Filtlong). DNA modification motif discovery and classification was performed with Nanodisco v1.0.3 (55), down-sampling read depth to 400× per plasmid, aligning reads to the pET15b vector containing the relevant mod insert. For genomic DNA libraries of E. coli and A. pleuropneumoniae, a 500 ng pool of 10 indexed samples was sequenced on 1 PromethION R9.4.1 flow cells (FLO-PRO002) on a PromethION 24 instrument running MinKNOW 22.08.6. Reads were basecalled in high-accuracy mode (dna_r9.4.1_450bps_hac_prom) with demultiplexing, using Guppy 6.1.5 and default parameters, and downsampled to a read depth of 400× per sample. DNA modification motif discovery and classification was performed with Nanodisco v1.0.3, aligning reads to the E. coli and A. pleuropneumoniae complete genomes. During the motif discovery step, a P-value threshold of –log₁₀(P) >20 was used for peak sequence selection. A threshold of –log₁₀(P) >100 was additionally used in analysis of A. pleuropneumoniae samples.

Overnight cultures of enriched ON/OFF variants (ModP1, ModP2 and ModQ) were grown to mid-log phase (OD₆₀₀ ∼ 0.4–0.7) and cultures were normalized to OD₆₀₀ of 1. Cells were harvested by centrifugation (5500 × g for 5 min) and resuspended in 300 μl of 6 M guanidinium chloride, 50 mM Tris–HCl (pH 8.0) and 10 mM dithiothreitol (DTT). Proteins were then alkylated by addition of 25 mM acrylamide and incubated for 60 min at 37°C with 500 rpm. Proteins were precipitated by addition of methanol:acetone (1:1) and kept overnight at –20°C. The protein was pelleted (18 000 × g, 10 min) and resuspended in 50 μl of 50 mM Tris–HCl containing 1 μg trypsin (NEB). After overnight incubation at 37°C, the trypsin digested peptides were cleaned with ZipTips (Millipore) and SWATH-MS analysis was performed as described previously (56). Briefly, tryptic peptides were analyzed by liquid chromatography-electrospray ionization-tandem mass spectrometry (LC–ESI-MS/MS) using a Prominence nanoLC system (Shimadzu) and a TripleTOF 5600 mass spectrometer with a NanoSpray III interface (Sciex). Peptides were separated on a Vydac Everest reversed-phase C18 high-performance liquid chromatography (HPLC) column at a flow rate of 1ml/min. A gradient of 10–60% buffer B over 45 min, with buffer A (1% acetonitrile and 0.1% formic acid) and buffer B (80% acetonitrile and 0.1% formic acid) was used. A mass spectrometry (MS)-time of flight (TOF) scan was performed from an m/z range of 350 to 1800 for 0.5 s, followed by information-dependent acquisition of MS/MS of the top 20 peptides from m/z 40 to 1800 for 0.05 s per spectrum, with automated CE (capillary electrophoresis) selection. Identical LC conditions were used for SWATH-MS. SWATH-MS of triplicate biological replicates was performed with the same MS-TOF scan, followed by high-sensitivity information-independent acquisition with m/z isolation windows with 1 m/z window overlap each for 0.1 s across an m/z range of 400–1250. Collision energy was automatically assigned by Analyst software (AB Sciex) based on m/z window ranges. Proteins were identified by searching against A. pleuropneumoniae strain AP76 and JL03 (NCBI accession no. CP001091.1 and CP000687.1, respectively) and common contaminants with standard settings using ProteinPilot 5.0.1 (AB Sciex). False-discovery-rate analysis was performed on all searches. ProteinPilot search results were used as ion libraries for SWATH analyses. The abundance of proteins was measured automatically using PeakView (AB Sciex) with standard settings. Comparison of protein relative abundance was performed based on protein intensities or ion intensities using a linear mixed-effects model with the MSstats package in R. Proteins with ≥2.0-fold changes in abundance and with adjusted P values greater than ≤0.05 were considered differentially expressed.

Enriched ON/OFF variants of A. pleuropneumoniae were grown to mid log phase (OD₆₀₀ ∼ 0.4–0.7) in BHI broth containing 0.01% NAD. The cells were centrifuged at 3000 × g for 15 min and resuspended in 10 mM Tris–HCl (pH 8.0). OMPs were prepared as described previously (57). Briefly, cells were sonicated and centrifuged at 6500 × g for 15 min. The supernatant was transferred to a new tube and sarkosyl was added to a final concentration of 1% (v/v). The mixture was incubated for 30 min at room temperature and ultracentrifuged at 15 000 × g for 90 min. The pellets were resuspended in 10 mM Tris–HCl (pH 8.0) containing sarkosyl (1% v/v final concentration) and incubated for 30 min at room temperature and, then ultracentrifuged at 15 000 × g for 90 min. The supernatant was discarded and the OMP-enriched pellet was resuspended in 100 μl of 10 mM Tris–HCl (pH 8.0). The protein concentration was measured using a Pierce™ BCA protein assay kit (Thermo Scientific). A 5 μg aliquot of each OMP preparation was run using the Novex Bis–Tris pre-cast gel system with 1× MOPS running buffer. Ammoniacal silver staining was carried out to visualize gross differences of protein expression due to ModP and ModQ phase variation (58). Differences in protein expression of selected proteins were quantified with ImageJ software by comparing the selected band intensity to a control band in the respective ON vs OFF lanes.

The minimum inhibitory concentration (MIC) assay was performed by a standardised broth microdilution method in Muller Hinton broth (MHB) with 2% yeast extract and NAD (0.1%) (59). The antimicrobials used were ampicillin, penicillin, florfenicol, ceftiofur, tulathromycin, tilmicosin and tiamulin. A. pleuropneumoniae ATCC 27090 was used as the quality control strain as recommended by CLSI guidelines (59). Briefly, overnight cultures of A. pleuropneumoniae enriched mod ON-OFF strains were grown to mid-log phase (OD₆₀₀ ∼ 0.4–0.7). Cultures were then diluted to an OD₆₀₀ of 0.2 and 50 μl of cultures were loaded into 96-well plates containing two-fold serially diluted antibiotics (ampicillin, penicillin, florfenicol from 5 to 0.1 μg/ml; tiamulin, tulathromycin, tilmicosin from 128 to 0.25 μg/ml; and enrofloxacin from 2 to 0.003 μg/ml). Plates were incubated for 24 h at 37°C under 5% CO₂. The MIC (μg/ml) was determined as the lowest concentration of antibiotic to inhibit bacterial growth. Each assay was performed in biological triplicate.

A growth curve assay was performed by diluting overnight cultures of enriched ON/OFF variants (modP1, modP2 and modQ) in BHI-NAD broth to OD₆₀₀ of 0.2. The cultures were then inoculated in triplicate into a 96-well plate and incubated at 37°C with shaking. Absorbance was measured (OD₆₀₀) every 15 min using a BMG Labtech Fluostar Optima plate reader for 18 h. Biofilm formation assays was performed as described previously (60). Enriched ON/OFF variants of modP1, modP2 and modQ were grown overnight on a BHI-NAD agar plate (0.01% NAD) at 37°C. The bacteria were harvested and suspended in fresh BHI-NAD broth and standardized to an OD₆₀₀ of 0.2, and 100 μl of this suspension was inoculated in triplicate into a 96-well plate. After overnight incubation at 37°C, the broth was removed from each well and the well washed with distilled water. A 100 μl aliquot of crystal violet (0.1%) was then added into each well and incubated for 15 min at room temperature. The wells were washed with water three times and dried for 30 min at 37°C. Biofilms were disrupted with 100 μl of 70% ethanol and absorbance was measured at OD₅₉₀ using a BMG Labtech Fluostar Optima plate reader.

GraphPad Prism 5.0 (GraphPad Software, La Jolla, California) was used to generate graphs and statistics. Error bars represent standard deviation from mean values. P values were considered as significant at <0.05 (*), <0.01 (**) and <0.001 (***), assessed using Student's t test (unpaired, two tailed). P values of <0.05 were considered not statistically different (NSD).

Our analysis of REBASE (22) revealed the presence of multiple Type I and Type III methyltransferases in A. pleuropneumoniae genomes with the hallmarks of phase-variation (23,25): duplicated variable hsdS genes containing IRs, or simple DNA sequence repeat (SSR) tracts. We observed the presence of multiple, duplicated hsdS genes in a Type I R-M system, with these duplicated variable hsdS genes containing inverted repeats (IRs). This Type I system also contained a gene annotated as a recombinase/integrase. The hsdR and hsdM genes encoded by this system were highly conserved (>95% nucleotide identity) in all 15 strains in REBASE identified to contain this system (25). The observed Type I system may therefore phase-vary much like the well characterised phase-variable Type I R-M system in S. pneumoniae strain D39 (34), the SpnD39III system, and subsequently described in multiple additional strains of S. pneumoniae (38,61). The prototype A. pleuropneumoniae strains AP76 (NCBI accession no. CP001091.1, locus tags APP7_1526-APP7_1531) and JL03 (NCBI accession no. CP000687.1, locus tags APJL_1433- APJL_1438) both encoded this system.

All the Type III mod genes we identified contained a simple DNA sequence repeat tract within the open reading frame of each mod gene (23). All previously characterised mod genes containing SSRs have been shown to be phase-variable, and to control a phasevarion (26,28–30,40). Analysis of the sequences of the Type III mod genes we identified in A. pleuropneumoniae demonstrated the presence of two separate, unrelated phase-variable mod genes, which we named modP and modQ. Both modP and modQ genes contained a GCACA_(n) SSR tract in their 5′ region, downstream of the annotated ATG start codon. Further analysis of these sequences demonstrated that modQ was represented by a single variant. The modQ TRD was the same in all strains of A. pleuropneumoniae in GenBank with a closed, annotated genome (28 strains in total) where a modQ gene is present (e.g. in strain JL03, NCBI accession no. CP000687, locus tag APJL_0819). Analysis of all modP sequences from these same 28 strains revealed the presence of four distinct TRD sequences, indicating at least four separate allelic variants of the modP gene present in the A. pleuropneumoniae population, similar to that seen with multiple other mod genes (26,28,39,62). We named these alleles modP1 (e.g. in strain AP76, NCBI accession no. CP001091.1, locus tag APP7_0747), modP2 (e.g. in strain JL03, NCBI accession no. CP000687.1; locus tag APJL_0704), modP3 (e.g. in strain S1536, NCBI accession no. CP031875.1, locus tag D1095_03465), and modP4 (e.g. in strain Femø, NCBI accession no. CP069796.1, locus tag D1099_04110).

Using conserved sequences of the identified Type I R-M system (the hsdR and hsdM regions), and each individual phase-variable Type III mod gene, we carried out a detailed sequence analysis using a diverse collection of 182 A. pleuropneumoniae isolates, and the 28 strains with a publicly available fully annotated whole genome to determine the distribution of each of these methyltransferases in A. pleuropneumoniae (Figure 2) (63). All 210 strains used for phylogenetic analysis are listed in Supplementary Table S1. This analysis demonstrated that all strains contained at least one of the Type III mod genes containing an SSR tract (19). Of these 210 strains, the modP1 allele was present in 150 strains followed by modQ in 61 strains, modP3 in 20 strains, modP2 and modP4 in 17 strains. Our analysis demonstrated that approximately 26% (n = 56) of strains harboured both modP and modQ genes. This analysis also revealed the distribution of Type III mod genes is serovar specific. For example, 100% of serovar 8 (n = 104) and serovar 7 (n = 21) encoded the modP1 gene whilst the modP3 gene was only present in serovar 2 (n = 20). All serovar 5 isolates (n = 6) were found to harbour the modP2 gene. For serovar 12, six isolates encoded the modP2 allele and one isolate encoded the modP4 allele, which was also present in all serovar 6 isolates (n = 10). The modQ gene appears to be absent in all isolates of serovars 4, 7, 8, 9, 17 and K2:07, as well as 3/21 serovar 2 isolates, whereas it is present in all other serovar isolates (Figure 2, Supplementary Table S1). In addition, 195 strains (93%) encoded the Type I locus encoding duplicate, variable hsdS loci containing inverted repeats (64). Of these 195 strains, 78% (n = 152) encoded both the Type I system and one of the Type III mod genes containing SSRs; 22% (n = 43) had both the phase-variable Type I system and two Type III mod genes containing SSRs. This analysis indicated that individual strains can contain both phase-variable Type I and Type III methyltransferases, a phenomenon never before observed in individual bacterial strains.

We also investigated the distribution of these newly identified Type III mod genes in a second collection of A. pleuropneumoniae strains, comprising 250 field isolates and 15 reference strains that represented one of each of the 15 serovars of A. pleuropneumoniae (65) (Supplementary Table S2). These field isolates comprise the Australian national reference collection for A. pleuropneumoniae (48,49). These field isolates have not been whole genome sequenced, and represented isolates from diverse geographical locations around Australia and across the globe. To determine the Type III genes/alleles present in this collection, we designed primer pairs specific for the TRDs of each modP allele, and for modQ, in order to screen via PCR. This analysis showed that the most prevalent serovars in Australia (serovars 1, 5, 7, 12 and 15) (66) harboured the same allelic variant(s) of these newly identified mod genes as in our first collection. For example, 100% of serovar 1 and serovar 15 isolates encoded the modP1 allele and the modQ gene. The modP1 allele and modQ gene were also present in A. pleuropnuemoniae reference strains of serovar 9 and 11. The modP2 allele and the modQ gene were found in all serovar 5 isolates (100%). For serovar 7, almost all (96%, n = 67) isolates were found to contain only the modP1 allele, with the exception of a single serovar 7 reference strain (WF83) and two Australian field isolates that contained the modP1 allele and modQ gene. The modP3 allele was present in serovar 2 (n = 2), 12 (n = 4), serovar 14 reference strain (n = 1) and two non-typeable strains. No strain in our culture collection encoded the modP4 allele except our single serovar 6 reference strain, Femø. Analysis of the field isolates revealed that 97% isolates encoded the modP gene (n = 256), with the modP1 allele being the most prevalent, present in ∼78% (n = 206) of all isolates. The modP2 allele was found at a lower frequency, comprising only 15% of isolates (n = 39) and only 3.8% (n = 10) isolates encoded the modP3 allele. In addition, 70% of A. pleuropneumoniae isolates in this collection (n = 185) contained the modQ gene. Overall, approximately 65% (n = 173) of isolates contained at least one of the four modP alleles and the modQ gene. Due to high prevalence of modP1, modP2 and modQ alleles in our Australian culture collections, we selected these three alleles for further analysis in this study.

In order to determine the extent of variable HsdS proteins that are encoded by the almost ubiquitous phase-variable Type I R-M system we have observed in A. pleuropneumoniae (Figure 2; illustrated in Figure 3A), we analysed the sequences of all TRD regions available in REBASE (15 strains). Analysis of these TRD sequences revealed that there are four unique 5′-TRD sequences, and four unique 3′-TRD sequences present in these strains (Figure 3B; full sequences and analysis in Supplementary Data 1), meaning a potential 16 unique HsdS proteins are possible dependent on the encoded TRDs present (four possible 5′-TRD x four possible 3′-TRD). This is unlike other phase-variable Type I systems, such as the SpnD39III system in S. pneumoniae (34,37,38,61), and the phase-variable Type I system in S. suis, (39) where all strains encode the same TRDs, and therefore shuffle between the same complement of expressed HsdS proteins. Our prototype strains of AP76 and JL03 encoded a subset of these sixteen potential TRDs: AP76 encoded a combination of 5’-TRD-1 and 3’-TRD-1 (abbreviated TRD1-1) and TRDs 1-2, 2-2 and 2-1 (comprising alleles A-D; Figure 3C); while in strain JL03 we found TRDs 2-1, 2-3, 3-3 and 3-1 (comprising alleles D-G; Figure 3D). In order to demonstrate that these hsdS alleles are phase-variable, we performed a semi-quantitative RT-PCR to demonstrate that all four potential hsdS allelic variants were expressed in a population of A. pleuropneumoniae, as carried out previously for a phase-variable Type I system in S. suis (39). This analysis demonstrated that all four hsdS alleles that we predicted to be present (alleles A-D in strain AP76, and alleles D-G in strain JL03) were expressed in a population of the respective strains (Figure 3C and D). This analysis provides excellent evidence that phase variation of this Type I system is occurring in both strains by shuffling of variable TRDs between the expressed and silent hsdS loci present. We also designed a FAM labelled PCR coupled to fragment length analysis to quantify the proportion of bacteria expressing each allele within a population of A. pleuropneumoniae (Supplementary Figure S1). Using these same prototype strains AP76 and JL03, we showed that the majority of the population in strain AP76 expressed hsdS allele A (5′-TRD-1 and 3′-TRD-1), and in strain JL03 hsdS allele E (5′-TRD-2 and 3′-TRD-3), but all other minor variants were also detectable using this methodology (Supplementary Figure S1D). Along with our semi-quantitative RT-PCR, this conclusively demonstrates phase variation of the Type I system in A. pleuropneumoniae.

As described above, we determined the presence of four unique modP alleles in A. pleuropneumoniae. Each allelic variant shows high sequence variability in the central TRD encoding region (illustrated in Figure 4A), with <25% nucleotide identity between each modP allele (Figure 4B), consistent with that seen for other phase-variable mod genes (39). The 5′ and 3′ regions of these modP alleles share >95% nucleotide sequence identity.

To examine if these modP genes were phase-variably expressed, we performed single colony isolation and enrichment to isolate individual variants of prototype strains encoding the two most common alleles, modP1 and modP2. This would allow enrichment for a triplet set of isogenic strains that were ∼90% pure for a single GCACA_(n) SSR tract length in the respective modP gene. We performed a FAM-labelled PCR coupled to fragment length analysis over the GCACA_(n) SSR tract to isolate these variants, each representing one of the three reading frames possible (Figure 4C). We used strain AP76 as our prototype strain for modP1 and strain JL03 for modP2. This resulted in a triplet of modP1 strains in AP76 where the GCACA_(n) SSR tract was 10, 11, or 12 repeats and for modP2 in JL03 with 23, 24 or 25 repeats (Figure 4C). Using whole cell lysates of these enriched triplet sets of each prototype strain we carried out Western blotting with antisera raised against the conserved region of ModP. This demonstrated that each ModP variant is only expressed when the GCACA_(n) SSR tract maintains the open reading frame (ON), and not expressed when the SSR tract length results in a frameshift and premature stop codon (OFF) (Figure 4D, full blot in Supplementary Figure S2). ModP1 in strain AP76 was only expressed when there were 10 GCACA repeats present, and ModP2 in strain JL03 when there were 25 GCACA repeats present.

Using long read sequencing from both PacBio SMRT sequencing and Oxford Nanopore technology, we performed methylome analysis of both ModP1 and ModP2. We did this by sequencing both plasmid vectors and genomic DNA isolated from E. coli BL21 strains over-expressing each cloned methyltransferase, and genomic DNA isolated from each of our triplet sets of prototype A. pleuropneumoniae strains encoding either the modP1 gene (strain AP76) or the modP2 gene (strain JL03) and modP1 and modP2 knockout A. pleuropneumoniae mutants. This analysis demonstrated that both ModP1 and ModP2 are cytosine specific Type III DNA methyltransferases, each recognising and methylating a different DNA target sequence. ModP1 methylates a cytosine in the motif CAA^(m4)CT, and ModP2 methylates a cytosine in the motif GWC^(m4)CT. Cytosine specific Type III DNA methyltransferases have been previously identified (67), but the ModP methyltransferase is the first identified phase-variable cytosine-specific Type III DNA methyltransferase (Table 2; all methylome data is presented in Supplementary Data 2; examples of PacBio and Nanopore analysis for each type of DNA analysed and motifs detected is presented in Supplementary Figures S3-S6).

Following this novel finding, we wanted to determine if ModP1 and ModP2 phase-variable expression resulted in the regulation of different phasevarions, i.e. did cytosine methylation at different motifs result in different sets of proteins being regulated commensurate with the observed phase variation? To quantify the extent of differential protein expression mediated by ModP1 and ModP2 phase variation, we carried out quantitative SWATH proteomics (68), used previously to characterise expression differences mediated by phase-variable DNA methyltransferases (40). The expression profiles of the ModP1 ON-OFF pair (10 repeats ON and 11 repeats OFF) and ModP2 ON–OFF pair (25 repeats ON and 24 repeats OFF) were assessed by SWATH-MS proteomics. Both strain pairs showed unique differences in protein expression dependent on which ModP protein was expressed, i.e. ModP1 and ModP2 control different phasevarions. Our SWATH-MS covered 26% of total proteins in strain AP76 (ModP1) and 32% of total proteins in strain JL03 (ModP2). In addition, gross protein expression differences were observed using outer membrane protein (OMP) fractions prepared from our enriched modP1 and modP2 ON-OFF pairs (Supplementary Figure S7).

In strain AP76 encoding ModP1, a total of 42 proteins had an expression ratio of ≥2.0 fold different between the ModP1 ON and ModP1 OFF isogenic strains. In ModP1 ON, 27 proteins were upregulated, and 15 proteins were downregulated (Table 3). Of the proteins showing increased expression in our ModP1 variant, many are involved in pathobiology, including the major RTX toxin, ApxIVA, a heat shock-like protein 33 HslO, and outer membrane protein W. Several proteins involved in energy metabolism were also upregulated in ModP1 ON including succinyl-CoA synthetase beta chain and mannose permease IID. In addition, several regulatory proteins were upregulated in ModP1 ON, such as leucine-responsive regulatory protein (Lrp) (69), and a putative HTH-type transcriptional. Several biosynthetic proteins showed an increase in expression in ModP1 ON, including those involved in amino acid biosynthesis (ketol-acid reductoisomerase, D-3-phosphoglycerate dehydrogenase, aspartate-semialdehyde dehydrogenase), heme biosynthesis (coproporphyrinogen III oxidase) and lipid A biosynthesis (acylUDP-N- acetylglucosamine O-acyltransferase). In the ModP1 OFF strain, many ribosomal proteins, outer membrane protein precursor PalA, lipoprotein and the fatty acid biosynthesis protein (malonyl CoA-acyl carrier protein transacylase (MCT)) exhibited increased expression compared to ModP1 ON (Table 3).

In strain JL03, phase variation of ModP2 resulted in varied expression of a different set of proteins, with 35 proteins showing a ≥2.0-fold expression ratio between ModP2 ON and ModP2 OFF, of which 9 proteins were upregulated in ModP2 ON, and 26 proteins downregulated in ModP2 ON (Table 4). Our ModP2 ON strain showing increased expression of proteins involved in central metabolism, including an acyl carrier protein, a biotin carboxyl carrier protein of acetyl-CoA carboxylase, a GTP binding protein and a transcription factor, NusA. In ModP2 OFF, proteins showing increased expression included the major RTX toxin, ApxIIIA, as well as a different set of proteins involved in general metabolism compared to ModP2 ON, such as glucose-6-phosphate 1-dehydrogenase, 2-dehydro-3-deoxyphosphogluconate aldolase and several amino acid biosynthesis proteins. The expression of several CRISPR associated proteins, three stress response factors (a large-conductance mechanosensitive channel, an excinuclease ABC subunit A, and sigma-factor 32) were also upregulated in ModP2 OFF, as were several proteins involved in fatty acid and phospholipid metabolism (putative acyl CoA thioester hydrolase, long-chain-fatty-acid–CoA ligase). ModP2 OFF also resulted in increased expression of key outer-membrane components, including a galactosyltransferase involved in lipooligosaccharide biosynthesis.

We also characterised the second observed mod gene, modQ that contained a GCACA_(n) repeat tract in A. pleuropneumoniae, (23) (Figure 5A). In order to determine if the modQ gene was also phase-variable, we also isolated a triplet set of strains with three consecutive repeat tract lengths in strain JL03 (GCACA_(n) repeat tracts of 9, 10, and 11 repeats; Figure 5B). Using ModQ specific antisera, we demonstrated that ModQ was only expressed when there were 9 GCACA repeats present in the SSR tract (Figure 5C, full blot in Supplementary Figure S8), and was OFF (not expressed) when either 10 or 11 GCACA repeats were present (Figure 5C). As our prototype strain JL03 encoded both ModQ and the ModP2 allele, we carried out Western blotting with our ModP2 and ModQ triplet set of strains to ensure that the mod gene under study contained the second mod gene in the same ON-OFF state, i.e, did our modQ ON-OFF pair contain the modP2 gene in the same state, and vice-versa. These Western blots (Supplementary Figure S9) demonstrated that in our ModQ ON-OFF triplet, where ModQ is either 9 GCACA repeats ON, or 10 or 11 GCACA repeats OFF (anti-ModQ antisera), modP2 was OFF in all three strains (anti-ModP sera with the same ModQ cell lysates). Analysis of the ModP2 triplet strain set showed that ModQ was ON in the ModP2 strain where there were 23 GCACA repeats, but ModQ was OFF in strains encoding ModP2 25 repeats ON, or ModP2 24 repeats OFF (anti-ModQ sera with the same ModP2 cell lysates). Therefore, for phenotypic analysis, we would use ModP2 ON GCACA 25 repeats, and ModP2 OFF GCACA 24 repeats (as described above).

Methylome analysis using genomic DNA from our triplet set of A. pleuropneumoniae JL03 strains encoding ModQ ON or OFF, and plasmid and genomic DNA from over-expression of ModQ and modQ knockout A. pleuropneumoniae mutant strain, demonstrated that the ModQ methyltransferase was an adenine specific Type III DNA methyltransferase, recognising the motif AG^(m6)ATG (Table 2; Supplementary Data 2). This was confirmed using both PacBio SMRT sequencing and Oxford Nanopore long read sequencing and methylome analysis (Table 2; Supplementary Data 2; Supplementary Figures S3–S6).

Quantitative SWATH proteomics with a ModQ ON-OFF pair (9 repeats ON and 10 repeats OFF) demonstrated a total of 86 proteins had an expression level of ≥2.0-fold difference between ModQ ON-OFF pair, with 62 proteins upregulated and 24 proteins downregulated in expression in ModQ ON strain (Table 5). In ModQ ON strain, several proteins were upregulated including major RTX toxin, ApxIIIC, outer membrane lipoprotein VacJ, cell envelop proteins and several stress response proteins such as excinuclease ABC subunit A (UvrA), conserved glutaredoxin-like protein, peptide methionine sulfoxide reductase, thioredoxin (TrxA) and survival protein E (SurE). Proteins involved in metabolic pathway as well as amino acid biosynthesis proteins were also showed higher expression in ModQ ON strain. A DNA biosynthetic enzyme, as well as proteins involved in nutrient uptake and energy metabolism showed increased expression in ModQ ON strain compared to ModQ OFF strain. In the ModQ OFF strain, several ribosomal proteins, and proteins involved in fatty acid and phospholipid metabolism (probable biotin carboxyl carrier protein of acetyl-CoA carboxylase, accB) were all upregulated. Three major outer membrane proteins were also upregulated in ModQ OFF strain compared to ModQ ON strain (OmpH family outer membrane protein, MomP1 and ModP2). Two cell division proteins (ZipA and FtsZ) also showed increased expression in ModQ OFF compared to ModQ ON strain (Table 5). In addition, examination of OMP fractions from ModQ ON and ModQ OFF showed gross differences in outer-membrane protein expression (Supplementary Figure S7).

We conducted standard growth curves of our enriched mod ON/OFF variants (modP1, modP2 and modQ) in BHI-NAD broth to determine if switching of mod genes could affect growth rates of A. pleuropneumoniae. Strain AP76 expressing ModP1 (ModP1 ON) exhibited a significantly higher growth rate compared to its isogenic strain that did not express modP1 (ModP1 OFF) (Figure 6A). Phase-variable expression of both modP2 and modQ had no impact on growth rate (Figure 6B and C). A standard static biofilm assay found that biofilms formed by the modP1 ON variant were larger in mass than its isogenic modP1 OFF variant (Figure 6A; P ≤ 0.01 determined with Student's t-test). No differences in biofilm formation were found between either the ModP2 or the ModQ ON versus OFF variants (Figure 6B and C), indicating that both growth rate and biofilm formation are influenced by ModP1 phase-variation.

Some human-adapted pathogens showed differential susceptibility to antibiotics due to phase variation of mod genes (28,70). Therefore, we aimed to demonstrate if any of the three phasevarions (ModP1, ModP2 or ModQ) influenced susceptibility to antibiotics regularly used to treat A. pleuropneumoniae infections in pigs. MIC is defined as the lowest concentration of antibiotic used to inhibit bacterial growth. MIC assays of our enriched ON/OFF pairs of three mod genes (modP1, modP2 and modQ) were carried out to investigate susceptibility to ampicillin, penicillin, florfenicol, ceftiofur, tulathromycin, tilmicosin and tiamulin. This analysis showed that whilst the ModP1 and ModQ phasevarions did not result in any differential responses to the antibiotics tested, expression of ModP2 (ON) resulted in significantly higher resistance (higher MIC) to all tested antibiotics except enrofloxacin compared to ModP2 OFF (Table 6). These findings suggested that global DNA methylation by ModP2 phasevarion has an impact on antibiotic susceptibility by as yet unidentified mechanisms.