M. bovis strain HB0801 (GenBank accession no. NC_018077.1) was isolated from the lung of infected beef cattle in Hubei Province, China, and its genome was fully sequenced by our laboratory (Qi et al., 2012). The M. bovis HB0801 attenuated strains HB0801-P115, HB0801-P150, and HB0801-P180, abbreviated as P115, P150, and P180, respectively, which exhibit progressively decreasing virulence, were used (Zhang et al., 2014). All the strains were propagated in pleuropneumonia-like organism (PPLO) medium supplemented with 10% horse serum (Thermo Fisher Scientific, Waltham, MA, USA) at 37°C for 48–72 h as described previously (Zhang et al., 2014).

The DNA of strains P115, P150, and P180 was extracted using bacterial genomic DNA extraction kits (Tiangen, Beijing, China). The 454 pyrosequencing method was used to determine the whole genome sequences of strains P115, P150, and P180. Three paired-end sequencing libraries with 8-kb inserts were constructed at the China Tianjin Biochip Corporation (Tianjin, China). For each sample, one-fourth of a PicoTiterPlate was run on a Roche/454 GS FLX sequencer (454 Life Sciences, Branford, CT, USA) using titanium chemistry according to the manufacturer's recommendations. Finally, 83,750,973 bases with 286,408 reads were obtained for strain P115, while 80,961,239 bases with 277,020 reads were obtained for strain P150, and 83,589,289 bases with 260,056 reads were obtained for strain P180, resulting in 85.7-fold (P115), 82.8-fold (P150), and 85.5-fold (P180) depths of sequencing. All the reads for each genome were assembled de novo by the GS De Novo Assembler (version 2.6). Approximately 95% of reads were assembled for each genome, resulting in 11 scaffolds with 52 non-redundant contigs for strain P115, three scaffolds with 49 non-redundant contigs for strain P150, and one scaffold with 53 non-redundant contigs for strain P180. All the scaffolds were ordered and oriented according to the genome architecture of strain HB0801 (Qi et al., 2012). The N50 contig lengths of the large contigs (>1 kb) of each strain were 31,437 bp (P115), 31,456 bp (P150), and 31,042 bp (P180). The total numbers of base pairs of the non-redundant contigs were 916,922 bp (P115), 919,584 bp (P150), and 918,470 bp (P180), which agrees with the previously reported genome size of 991,702 bp for strain HB0801. To fill gaps within the scaffolds, polymerase chain reactions (PCRs) and capillary electrophoresis sequencing were performed with primers designed near the gaps. After gap filling, all the generated reads were respectively mapped to the corresponding full genome sequence using Burrows–Wheeler Aligner software (Li and Durbin, 2010). With the help of deep sequencing coverage, possible homopolymer errors resulting from the 454 sequencing method were rectified manually.

Open reading frames of the three passaged genomes (P115, P150, and P180) were predicted initially using Glimmer 3.02 (http://www.cbcb.umd.edu/software/glimmer/) and then modified using the translated nucleotide Basic Local Alignment Search Tool (BLAST) algorithm (http://blast.ncbi.nlm.nih.gov/) and compared to the genome of the primary strain HB0801 (GenBank accession no. CP002058). The functions of the coding sequences in the passaged strains were defined by referencing strain HB0801 (Qi et al., 2012). A comparative analysis between these four strains (P115, P150, P180, and HB0801) was conducted using MEGA 5.0 (Tamura et al., 2011) and Mauve 2.3.1 (Darling et al., 2004) genome alignment software. In addition, we detected differences in SNPs and indel sites among these strains.

To confirm the presence of the 14.2-kb deleted region during continuous passage, a putative lipoprotein-encoding gene (Mbov_0732) present in the deleted region was selected for PCR with forward (5′–AGCGACCAAAATACTAGAC –3′) and reverse (5′–TCGTTGCCACTGTATTCA–3′) primers using the following program: 95°C 3 min, followed by 30 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C 2 min, followed by 72°C for 15 min and 16°C for 5 min.

SNPs and indels might affect the expression of essential genes. Sixty-seven pairs of primers specific to the flanking sequences of the SNPs and indel sites were designed (Table S1) and Sanger DNA sequencing was performed with an ABI 3730 sequencer by China Tianjin Biochip Corporation (Tianjin, China).

Differentially expressed proteins predicted by genome comparisons of the virulent HB0801 strain to its attenuated P115, P150, and P180 derivatives were classified into two categories: (i) proteins deleted in P115, P150, or P180 (Table 1) and (ii) proteins displaying non-sense SNPs and indels (Table 2). Protein functions and related metabolic pathways were assigned according to the Uniprot database (http://www.uniprot.org/) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.kegg.jp/kegg-bin/show_organism?menu_type=pathway_maps&org=mbi), respectively. To determine whether proteins are secreted and contain a signal peptide, we used Pred-lipo (http://bioinformatics.biol.uoa.gr/PRED-LIPO/input.jsp), the SignalP 4.1 server for the prediction of classical secreted proteins (http://www.cbs.dtu.dk/services/SignalP/), and the SecretomeP 2.0 server for the prediction of non-classical secreted proteins (http://www.cbs.dtu.dk/services/SecretomeP/). Moreover, to assess whether the SNPs were present at active sites or domains of the proteins and to obtain the conserved domains in the secretory proteins, we used the National Center for Biotechnology Information (NCBI) conserved domain database (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). Paralogs of the genes and percentages of identity between paralogs were determined with the NCBI BLAST algorithm (http://blast.ncbi.nlm.nih.gov/Blast.cgi). Protein–Protein interactions were obtained using the STRING database (www.string-db.org).

The involvement of these proteins in virulence was further confirmed based on the whole genome of M. bovis HB0801 using the VFDB (http://www.mgc.ac.cn/VFs/) (Chen et al., 2016). The amino acid sequences of all the proteins were obtained from the NCBI. Each protein was aligned individually against the VFDB full dataset by the BLAST algorithm. A matrix was created by VFDB output, including VFDB hits against each protein in M. bovis, a related BLAST score, and an E-value. The matrix was filtered based on a BLAST score ≥80.

H₂O₂ production was measured as described previously in Mycoplasma mycoides subsp. capri (Allam et al., 2012). Briefly, eight transposon-disrupted mutants in the deleted region of the attenuated strains, including Mbov_0723, Mbov_0724, Mbov_0725, Mbov_0727, Mbov_0730, Mbov_0732, Mbov_0734, and Mbov_0735, as well as different passages of M. bovis HB0801, including P1, P25, P50, P75, P100, P115, P150, and P180, were grown to mid-logarithmic phase. The transposon-disrupted mutants were prepared using the wild-type HB0801 strain and a Tn4001 transposon. The transposon disrupted the abovementioned genes that are present in the deleted region mentioned in Table 1. Moreover, the Mbov_0723, Mbov_0724, Mbov_0725, Mbov_0727, Mbov_0730, Mbov_0732, Mbov_0734, and Mbov_0735 genes were transposon-disrupted at positions 85,6071, 85,7345, 858,521, 860,287, 863,914, 864,069, 866,042, and 866,786, respectively, in the genome. To assess H₂O₂ production, bacteria were harvested by centrifugation at 15,400 g for 4 min at 4°C. The pellets were washed three times in incubation buffer (67.7 mM HEPES, pH 7.3; 140 mM NaCl, and 7 mM MgCl₂). After the final wash, the M. bovis pellets were resuspended in incubation buffer to a density of 10⁸ cells/ml and incubated for 20 min at 37°C. After the incubation, 10 μl of glycerol (10 mM) was added, and the samples were kept at 37°C for 30 min. H₂O₂ production was determined using commercial H₂O₂ assay kits (Cayman Chemical, Ann Arbor, MI, USA). Standard curve determination and specificity quality controls were performed according to the instructional manual. The experiments were performed in three independent replicates.

To determine the statistical significance among the differences in SNPs after comparing all the genes of M. bovis, as well as among different SNP and indel types, a z-test for proportions (z = pˆ-p0/√p0(1 − p0)/n) was applied to generate a p-value. Moreover, to determine the statistical significance of the genes involved in different KEGG pathways, a KEGG enrichment analysis was performed using the KOBAS 2.0 web server, and Fisher's exact test was used to analyze the data (Xie et al., 2011). To analyze differences in H₂O₂ production among the different strains, two-way analysis of variance was performed using GraphPad Prism version 5.0 (GraphPad Software, San Diego, CA, USA). p < 0.05 were considered to be significantly different and are marked with an asterisk in the figures, while p < 0.01 were considered to very significantly different and are marked with two asterisks in the figures.

The complete genome sequences of the three attenuated strains were deposited in the GenBank database with the accession numbers CP007589 for P115, CP007590 for P150, and CP007591 for P180.

Because of a lack of genetic tools and limited research techniques, there is limited understanding of the pathogenesis and virulence-related factors of M. bovis. Hence, the whole genomes of strains P115, P150, and P180 were sequenced and assembled, and a comparative study was performed to identify significant virulence-related factors. The genome structures of the virulent M. bovis HB0801 strain and the three attenuated strains are shown in Figure 1A. An analysis of the three attenuated strains and the wild-type HB0801 strain revealed that all the genomes of strains P115, P150, and P180 lost a 14.2-kb fragment compared with strain HB0801 (Figure 1B and Table 1). Further analysis of the deleted fragment showed that it consists of 14 putative genes encoding ascorbate-specific phosphotransferase system (PTS) system enzyme IIB and IIA components, a phosphotriesterase family protein, predicted hydrolases of the haloacid dehalogenase (HAD) superfamily, DNA methyltransferase, and a type III restriction-modification (RM) system methylase, as shown in Table 1. In addition, 38 strains that were passage between one and 115 times were shown by PCR to have a deletion in the Mbov_0732 gene. The results demonstrated that the 14.2-kb deletion occurred in passage 115 and onward as shown in Figure S1, indicating that the deletion was stably maintained in the passaged strains after its occurrence.

In addition to the 14.2-kb deleted region, a genome comparison between strain HB0801 and the attenuated strains revealed many SNPs and indels. Sixty-seven non-sense SNP sites were tested with PCR and DNA sequencing, and 38 SNPs were confirmed (Table 2). In agreement with the increasing attenuation tendency resulting from continuous passage, strains that were subjected to more passages had more SNPs, and the SNPs in the strains that were subjected to fewer passages were shared by the higher passaged strains. The numbers of SNPs were 31, 24, and 19 for the P180, P150, and P115 strains, respectively. In addition, eight nucleotide deletions were found, and their frequencies of occurrence followed a pattern that was similar to that of the SNPs (Table 2). These changes in the genomes were analyzed statistically with a z-test of proportions, and the results showed that there were significant differences in the numbers of SNPs, indels, and SNPs+indels between the wild-type strain and each of the attenuated strains, and between any two attenuated strains (normalized p < 0.00001).

Theoretically, the functions of some genes containing SNPs or deletions could be compensated by their paralogs in the genome. Hence, the FASTA sequences of such paralogs were aligned with a reference protein by the BLAST algorithm to determine the percentage of identity between them (Table S2); a reference protein is a protein that is found in our genomics data (Tables 1, 2). The results showed that the proteins that have paralogs in our genomics data included four of the 14 genes in the 14.2-kb deleted region (Mbov_0727, Mbov_0729, Mbov_0732, and Mbov_0856) and 14 of the 46 genes containing SNPs or indels (Mbov_0018, Mbov_0049, Mbov_0111, Mbov_0134, Mbov_0338, Mbov_0347, Mbov_0350, Mbov_0393, Mbov_0518, Mbov_0525, Mbov_0581, Mbov_0682, Mbov_0742, and Mbov_0832). Hence, the genes responsible for virulence attenuation would likely include the other 10 genes in the 14.2-kb deleted region and the 32 genes with SNPs and indels, but which lacked paralogs.

The proteins listed in Tables 1, 2 were input into the STRING database to predict potential protein-protein interactions using the type strain M. bovis PG45 as the database default reference. Protein–protein interaction maps are given for 13 of the 60 input proteins (Figure 2). Different colored lines show different types of interactions. Hence, more lines show more interactions and an increased probability of a protein-protein interaction. The interactions among the proteins are highlighted in evidence view (Figure 2A). In evidence view, all possible interactions are shown. Different colored lines show different types of interactions e.g., gene fusion, co-occurrence, co-expression, experiments, databases, text mining etc. While in interactive view, the interactive proteins are clustered together based on interaction and co-occurrence. The evidence (Figure 2A) and interactive views (Figure 2B) showed that the most significant proteins that are likely to interact include enolase (Mbov_0482), pyruvate kinase (Mbov_0155), transketolase (Mbov_0212), L-lactate dehydrogenases (Mbov_0565), and D-lactate dehydrogenases (Mbov_0160), suggesting that these proteins might function together, and thus, that they are the strongest candidates for explaining the attenuated virulence of highly passage M. bovis strains.

A KEGG enrichment analysis was performed to determine whether the mutated genes in Tables 1, 2 are enriched in some important metabolic pathways. Fifteen genes were shown to be enriched in 15 metabolic pathways. Furthermore, 10 of the 15 genes were involved in more than one pathway (Figure 3 and Table S3). The 10 genes and the involved numbers of pathways numbers are Mbov_0155 and Mbov_0482 (seven), Mbov_0206 and Mbov_0565 (six), Mbov_0212 (five), Mbov_0567 (four), Mbov_0338 (three), and Mbov_0522, Mbov_0722, and Mbov_0723 (two). Among them, Mbov_0155, Mbov_0212, Mbov_0482, and Mbov_0565 were also identified in the protein-protein interaction analysis. Theoretically, proteins involved in more pathways should be more likely to play a significant role in M. bovis virulence. However, among the pathways, biosynthesis of secondary metabolites, biosynthesis of antibiotics, carbon metabolism, pyruvate metabolism, biosynthesis of amino acids, glycolysis/gluconeogenesis, pentose phosphate pathway, purine metabolism, and ATP-binding cassette (ABC) transporters were enriched in more pathways (six, six, five, four, four, four, three, three, and three, respectively).

Twelve predicted secretory proteins with signal peptides, which includes those encoded by 10 genes (Mbov_0049, Mbov_0111, Mbov_0347, Mbov_0350, Mbov_0393, Mbov_0518, Mbov_0525, Mbov_0579, Mbov_0682, and Mbov_0797) that have SNPs (Table 2) and the Mbov_0729 and Mbov_0856 genes, which are deleted, were identified (Table 1). A conserved domain analysis (Table 4) identified the presence of conserved domains in the eight secretory proteins mentioned in Table 3, but not in Mbov_0111, Mbov_0393, Mbov_0525, and Mbov_0856. Liproproteins and secretory proteins might contribute to bacterial virulence and pathogenesis because they are involved in many phenomena, ranging from cellular physiology through immune responses to virulence. Some secretory proteins were discovered previously in M. bovis (Khan et al., 2016; Zhang et al., 2016). Many other genes in the M. bovis genome are predicted to encode putative lipoproteins and secretory proteins, but whether these proteins exist in a functional form and how they are related to virulence need to be determined. The most significant conserved domains with available functions include CDC45 (Mbov_0797), DUF1388 (Mbov_0797), structural maintenance of chromosomes (SMC) (Mbov_0049), Rad23 (Mbov_0682), and peptidase C19 (Mbov_0729) (Table 4). CDC45 is required for the initiation of DNA replication (Saha et al., 1998). Members of the DUF1388 family function as the main targets for neurofilament-directed protein kinases in vivo. SMC proteins are essential for successful chromosome transmission during replication, and segregation of the genome, including chromosome condensation, recombination, DNA repair, and epigenetic silencing of gene expression (Haering et al., 2002). Rad23 family proteins are used for targeting nucleotide excision repair to specific parts of the genome. Peptidase C19 is a deubiquitinating enzyme that can deconjugate ubiquitin or ubiquitin-like proteins from ubiquitin-conjugated proteins (De Jong et al., 2006). The functions of the other domains are currently unknown.

To further investigate the virulence-related factors that contribute to the attenuation of highly passaged M. bovis strains, all the proteins of M. bovis were analyzed using the VFDB. In the VFDB full dataset, all proteins related to known and predicted virulence-related factors are present. Seventy-two genes in the M. bovis genome were shown to encode virulence-related factors based on a BLAST score ≥80 (Table 5). Although no virulence-related factor-encoding genes were found in the 14.2-kb deleted region, eight genes overlapped with those containing SNPs and indels (Table 2). In decreasing order of BLAST scores, they are Mbov_0482, Mbov_0533, Mbov_0581, Mbov_0338, Mbov_0134, Mbov_0742, Mbov_0018, and Mbov_0797.

By combining the above results for the deleted genes and the genes with SNPs and indels, 11 critical genes that likely contribute to the attenuation of highly passaged M. bovis strains were identified, and they include (in order of decreasing importance) Mbov_0722, Mbov_0723, Mbov_0482, Mbov_0565, Mbov_0155, Mbov_0581, Mbov_0742, Mbov_0299, Mbov_0212, Mbov_0797, and Mbov_0567.

Among the deleted genes, Mbov_0722 and Mbov_0723 are highlighted. They encode ascorbate-specific PTS enzyme IIB and IIA components, respectively, which function in ascorbate and aldarate metabolism pathways by involving the PTS, the major carbohydrate transport system in bacteria, and they are responsible for the conversion of L-ascorbate into L-ascorbate-6 phosphate, as shown in Figure S2 (Postma et al., 1993). The final product of this pathway is D-xylulose-5 phosphate, which is an intermediate in the pentose phosphate pathway. As is known, the primary purpose of this pathway is to generate a reducing equivalent of NADPH that can be used in reductive biosynthesis reactions, while the production of the pentose phosphate pathway intermediates ribose 5-phosphate and erythrose 4-phosphate are used to synthesize nucleotides and nucleic acids, and aromatic amino acids, respectively. In addition, this sugar has a role in gene expression, mainly by promoting the ChREBP transcription factor in the well-fed state (Iizuka and Horikawa, 2008). Hence, because of the deficiencies of both proteins in the attenuated strains, these related metabolic functions would be less efficient. Because the occurrence of the 14.2-kb deleted region began at passage 115 and was maintained from passages 115 through 180, these metabolic defects might contribute to the increasing attenuation of highly passaged M. bovis strains.

The enolase encoded by Mbov_0482 is considered to be a virulence-related factor, and it participates in seven metabolic pathways and interacts with other proteins as described previously. It catalyzes the reversible conversion of 2-phosphoglycerate into phosphoenolpyruvate. This reaction is present in the glycolysis/gluconeogenesis pathways of M. bovis, as shown in Figure S3. Prokaryotic α-enolase may contribute to pathophysiological processes (Pancholi, 2001). A surface-associated enolase is an adhesion-related factor of M. bovis that contributes to adherence by binding plasminogen (Song et al., 2012). The immunogenicity of enolase has also been observed in Mycoplasma synoviae (Bercic et al., 2008) and Mycoplasma capricolum subsp. capripneumoniae (Zhao et al., 2012). Hence, enolase might be a significant protein that contributes to the virulence of M. bovis. However, two SNPs were not detected inside the enolase-encoding gene, but were located approximately 300 bp upstream of the gene, suggesting that both SNPs might cause some change in the promoter region.

Like enolase, L-lactate dehydrogenase encoded by Mbov_0565 is another protein that is important for metabolism and as a virulence-related factor. It interconverts L-lactate to pyruvate. Additionally, this enzyme works in many metabolic pathways, including glycolysis/gluconeogenesis, cysteine and methionine metabolism, pyruvate metabolism, and propanoate metabolism. Hence, this lactate dehydrogenase plays significant roles in a variety of metabolic processes. Moreover, L-lactate dehydrogenase was shown to be surface expressed and to interact with plasminogen (Grundel et al., 2015). Furthermore, because of the SNP in the L-lactate dehydrogenase-encoding gene, a polar amino acid (threonine) was converted to a non-polar amino acid (methionine).

Pyruvate kinase (Mbov_0155) is also a very active protein in metabolism and energy production. It catalyzes the conversion of phosphoenolpyruvate to pyruvate (Figure S4). Then, mycoplasmas generate ATP via a proton-translocating ATP synthase by oxidizing organic acids (pyruvate and lactate) to acetate and CO₂. Moreover, the pyruvate kinase present in the genome of Mycoplasma suis is proposed to be required for the conversion of all NDPs and dNDPs to NTPs and dNTPs, respectively (Pollack et al., 2002). Hence, this protein is very important for energy production in Mycoplasma species, and the SNP in the pyruvate kinase-encoding gene may lead to the attenuation of M. bovis.

The different molecular functions of transketolase (Mbov_0212) include metal ion binding and transferase activity, and it is a key enzyme in the non-oxidative branch of the pentose phosphate pathway that transfers a two-carbon glycolaldehyde unit from a ketose donor to aldose-acceptor sugars (Jores et al., 2009). Moreover, transketolase is an immunodominant membrane protein that may be used as a biomarker for the serological diagnosis of contagious agalactia caused by M. mycoides subsp. capri (Corona et al., 2013). Hence, this protein may be involved in the immunogenicity and virulence of M. bovis.

ABC transporter proteins have different functions, including ATP binding, ATPase activity, and catalyzing the transmembrane movement of substances, such as importing sugars, amino acids, peptides, metal ions, and phosphates, and effluxing toxins, drugs, and proteins (Higgins et al., 1986). The Mbov_0742 gene, which contained a SNP only in strain P180, encodes a glycerol ABC transporter protein. The Mbov_0581 gene, which had an indel only in strain P180, encodes a multiple sugar ABC transporter protein. Both genes are considered to encode virulence-related factors in the VFDB. In addition, the Mbov_0134 and Mbov_0018 genes encoding ABC proteins responsible for the uptake of simple sugar and spermidine/putrescine were found to have SNPs in all three attenuated strains.

NADH dehydrogenase (Mbov_0299) is considered to be a virulence-related factor. It generates energy by transferring electrons from NADH (oxidation) to quinine. NADH and NADPH are possibly essential for the growth of M. suis (Guimaraes et al., 2011). In Mycobacterium tuberculosis, a mutant lacking NuoG, a subunit of the type I NADH dehydrogenase complex, exhibited attenuated growth in vivo (Blomgran et al., 2012). Hence, this gene mutation might be significantly related to the further attenuation of M. bovis strain P180.

Phosphate acetyltransferase (PTA) (Mbov_0567) catalyzes the formation of acetyl phosphate and acetyl CoA, which are used in many metabolic pathways, e.g., in taurine and hypotaurine metabolism, pyruvate metabolism, propanoate metabolism, and methane metabolism (Figure S5). Acetyl-phosphate regulates various cellular processes, including cell division, outer membrane protein expression, osmoregulation, and biofilm development. Moreover, acetyl phosphate is required for the activation of the Rrp2-RpoN-RpoS pathway, which serves as a global signal in bacterial pathogenesis by activating virulent genes (Xu et al., 2010). Moreover, Salmonella enterica serovar Typhimurium (Kim et al., 2006) and Vibrio cholerae (Chiang and Mekalanos, 1998) PTA mutants were shown to exhibit impaired growth and attenuated virulence. In addition, because of a SNP in Mbov_0567, a proline residue was converted to threonine, resulting in polarity change that might be significant in virulence attenuation.

In addition, membrane proteins influence cell shape, cell division, motility, and adhesion to host cells, and they are thought to be integrally involved in the pathogenesis of mycoplasmas. As is known, adhesion and invasion are generally considered to be virulence-associated processes. M. bovis can adhere to and invade epithelial cells and immune cells. Membrane lipoproteins, such as Vsps, enolase, and Vpmax, play significant roles in the adhesion of M. bovis to host cells (Burki et al., 2015). Five genes encoding membrane proteins were shown to have SNPs, and four of them were predicted to encode secretory proteins with signal peptides, including Mbov_0393, Mbov_0525, Mbov_0579, and Mbov_0797. Among them, Mbov_0797 encodes a Vsp, and it was predicted to be a virulence-related factor, Mbov_0579 was predicted to encode the ADP-ribosyltransferase CDTa, which contains functional domains of the community-acquired respiratory distress syndrome toxin of M. pneumoniae (Kannan et al., 2014).

H₂O₂ is thought to be a significant virulence related factor in M. pneumoniae, M. mycoides, and Mycoplasma ovipneumoniae. Secondary metabolites are considered to play significant roles in the pathogenesis of some Mycoplasma species (Pilo et al., 2005; Hames et al., 2009). For example, H₂O₂ was demonstrated to be a major virulence-related factor that leads to cell death and lipid peroxidation in M. pneumoniae (Hames et al., 2009) and MmmSC (Pilo et al., 2005). Moreover, activation of glycerol utilization and overproduction of H₂O₂ occurred during intracellular infection with Mycoplasma gallisepticum (Matyushkina et al., 2016). Although an in vitro H₂O₂ assay for M. bovis field strains showed that variations in H₂O₂ production did not correlate with M. bovis virulence (Schott et al., 2014), in vitro passaging an M. bovis strain resulted in decreased levels of H₂O₂ production (Khan et al., 2005).

To confirm that H₂O₂ production was affected by the mutations in the attenuated strains, H₂O₂ production was tested in the wild-type strain HB0801 (P1), strains of various passages (25, 50, 75, 100, 115, 150, and 180), and eight mutants specific to the genes in the 14.2-kb deleted region (Mbov_0723, Mbov_0724, Mbov_0725, Mbov_0727, Mbov_0730, Mbov_0732, Mbov_0734, and Mbov_0735). The results showed a decreasing tendency of H₂O₂ production following increasing numbers of passages in M. bovis. The differences in H₂O₂ production between passage 1 and the serially passaged strains began to be statistically significant from passage 115 and onward (p < 0.01) (Figure 4A). Interestingly, other researchers measured the production of H₂O₂ by in vitro passaged strains including the 50th, 100th, and 200th passages of M. bovis, and they showed a similar decreasing tendency of H₂O₂ production with increasing numbers of passages (Khan et al., 2005). Random transposon mutagenesis was used recently to generate M. bovis mutants (Sharma et al., 2014). With the availability and further improvement of these techniques, it should be possible to obtain detailed information about the interactions of M. bovis with its host in the near future. To differentiate the effects resulting from the 14.2-kb deleted region and the SNPs in related genes on H₂O₂ production, transposon-disrupted mutants were further characterized. However, there was no significant difference in H₂O₂ production by the mutants compared with the P1 strain (Figure 4B). Therefore, these genes in the 14.2-kb deleted region might not affect H₂O₂ production, or a single gene could only contribute slightly to H₂O₂ production. In fact, other genes outside the 14.2-kb deleted region that have SNPs and indels might be associated with altered H₂O₂ production via the impairment of various pathways (Figure 5), such as ABC transporters (Figure S6), carbon metabolism (Figure S7), biosynthesis of amino acids (Figure S8), glycolysis/gluconeogenesis (Figure S9), and pyruvate metabolism (Figure S10).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.