Bovine M. tb strain 1458 was isolated in this lab. The M. tb strain H37Rv and M. bovis strain AF2122_97 were kindly offered by Prof. Chuanyou Li from Beijing Tuberculosis and Thoracic Tumor Research Institute. The strains were cultured to mid-log phase in 20 ml of Middle brook 7H9 medium (Becton Dickinson and Company, Franklin Lakes, NJ, USA) that was supplemented with 10% oleic acid, albumin, dextrose, and catalase medium (Becton Dickinson and Company) and 0.05% Tween 80 (Amresco Inc., Solon, OH, USA), with agitation at 37°C for 1–2 weeks in a biosafety level 3 facility at Huazhong Agricultural University (Wuhan, China). The culture of M. tb strain 1458 was centrifuged at 10,000 × g and 4°C for 10 min, and DNA was extracted using cetyl trimethylammonium bromide, and 2 μg of DNA was used for sequencing. The total RNA of all the three strains was extracted using FastRNA Blue Kit (MP Biomedicals, Shanghai, China), and was reverse-transcribed into cDNA using a HiScript II One Step RT-PCR Kit (Vazyme Biotech, Nanjing, China).

Complete genomic sequencing was conducted using a Roche GS FLX system (454 Life Sciences Corp, Branford, CT, USA; Margulies et al., 2005). A total of 267,232 reads amounting 128 Mb of raw data (average read length of 478 bp), were obtained, resulting in 29-fold genome coverage. Assembly was performed using the GS de novo Assembler software from SFF file with default parameters, which will remove reads shorter than 50 bp. A total of 113 contigs ranging from 500 to 203,016 bp were produced. The relationships among the contigs were determined by multiplex polymerase chain reaction (PCR) (Tettelin et al., 1999). Then, gaps were filled in by sequencing the PCR products using ABI 3730XL capillary sequencers (Applied Biosystems, Waltham, MA, USA). Phred, Phrap, and Consed software packages (http://www.phrap.org/phredphrapconsed.html) were used for final assembly and editing, and the low-quality regions (the bases with phred quality score <20) of the genome were re-sequenced. The final sequencing accuracy was 99.99%.

Putative coding sequences were identified by Glimmer 3.02 (Delcher et al., 1999), and the peptides shorter than 30 amino acids were eliminated. Insertion sequences were first detected using the IS Finder database (https://www-is.biotoul.fr/) by the default parameters. The tRNA genes were predicted by tRNAScan-SE (Lowe and Eddy, 1997), and rRNA genes were predicted by RNAmmer (Lagesen et al., 2007) using the default parameters. Functional annotation of the coding sequences was performed by searching against NCBI non-redundant protein database using the protein Basic Local Alignment Search Tool (BLASTP) (E-value is set to 1e-10) (Altschul et al., 1997) and the Conserved Domain Database (CDD) (Marchler-Bauer et al., 2007) by reversed positon specific BLAST with E-value < 1e-5. The pseudogenes were identified by BLASTP, with an amino acid alignment length <20% of the referenced proteins' amino acid lengths. Metabolic pathways were constructed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2004). Genome comparisons were performed using Mauve (Darling et al., 2004). The genome atlas was drawn using GenomeViz1.1 (Ghai et al., 2004). A new pan-genome analysis pipeline (Zhao et al., 2012) was used to identify the orthologs among the three genomes, with the MP method using the following settings: coverage >50% and identity >50%.

Orthologs of known Mycobacterium tuberculosis complex (MTBC) genomes were obtained from the National Center for Biotechnology Information database (https://www.ncbi.nlm.nih.gov/genome/). The phylogenetic position of the genome of M. tb strain 1458 was determined based on orthologous proteins. Concatenated protein sequences of orthologous Mycobacterium species proteins were first aligned using MAFFT (Katoh and Standley, 2013), and then the conserved alignment blocks were extracted by the MEGA6 program (Tamura et al., 2013). A maximum likelihood tree was built with PHYML (Guindon et al., 2010) with 1,000 bootstrap replications and the following parameters: “JTT” for the substitution model; “estimated” for the proportion of invariable sites; “estimated” for the gamma distribution parameters; “4” for the number of substitution categories; “yes” to the optimized tree topology; and “BIONJ” for the starting tree.

PCR was used to confirm IS6110 insertion in the upstream of the genes esxS, phoP, dnaN, and ctpD in M. tb strain 1458. Primer sequences and products' size were listed in Table 1, and the PCR conditions were as follows: 95°C for 5 min, followed by 35 cycles of 95°C for 60s, 62°C for 60s, and 72°C for 90s, and ended with a final extension at 72°C for 10 min.

The absolute quantitative real-time PCR (qPCR) was used to determine the copy number of four typical genes (esxS, phoP, dnaN, and ctpD). The qPCR standards for these genes were firstly prepared. The primer sequences and products' size were listed in Table 2, and the four gene fragments were amplified and sequences verified. Then, the four PCR products were cloned into the pMD18-T vector (TaKaRa, Dalian, China) to construct the four recombinant plasmids, pMD-phoP, pMD-esxS, pMD-dnaN, and pMD-ctpD. The concentrations of the four recombinant plasmids were measured with a NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA), and the copy number was calculated as the following formula: number of copies = [recombinant plasmid concentration (ng) × 6.022 × 10²³]/[recombinant plasmid length × 1 × 10⁹ × 650].

The absolute qPCR was performed by using AceQ qPCR SYBR® Green Master Mix (Vazyme Biotech, Nanjing, China). The volume of each reaction was 10 μl, including 1.5 μl of cDNA, 0.25 μl of each primer, 5 μl of mixture and double distilled water to a volume of 10 μL. The whole experiment was programmed in an ABI ViiA™ 7 Real-Time PCR System (Applied Biosystems, Waltham, MA, USA) as follows: 95°C for 10 min for the hot-start, followed by 40 cycles of 95°C for 30 s, 58°C for 30 s, and 72°C for 30 s, and a melting curve program. The fluorescence signal was collected at the end of each elongation step. Ten-fold serial dilutions of the four recombinant plasmids DNA from 10¹ to 10¹⁰ copies/μl were prepared to generate a standard curve for absolute qPCR. The results are expressed as the mean ± the standard deviation (SD) of the mean. Differences in the data between the groups were analyzed by ANOVA using the GraphPad Prism5 and p-values < 0.05 and < 0.01 were considered to be statistically significant (^*) and very significant (^**).

The genome of the M. tb 1458 strain contains a single circular chromosome of 4,402,033 bp with a GC content of 65.6% (GenBank accession no: CP013475; Figure 1). We identified 4,349 open reading frames (ORFs) in the genome, with an average length of 915 bp, which occupy 90.4% of the genome. Among these ORFs, 2,898 (66.6%) genes could be classified into Clusters of Orthologous Groups (COGs) families comprising 21 functional categories (Figure 2, Table S1).

To perform a comprehensive comparative genomic analysis, a total of 25 mycobacterial genomes in addition to that of M. tb strain 1458 were retrieved from the GenBank database, and a phylogenetic analysis was performed. As expected, all four Beijing strains fell within the same cluster, and the 1458 strain was most closely related to the Beijing family Chinese strain CCDC5079 (Figure 3).

In total, 18 kinds of complete IS elements were found in the M. tb 1458 genome using online tools (https://www-is.biotoul.fr/). Basically, the IS elements belong to the IS256, IS607, IS3, and IS521 families (Table S2). In addition, we identified 21 copies of IS6110 in strain 1458, compared with 16 copies in strain H37Rv. Among them, 11 IS6110s are intragenic insertions (Table 3), and the disrupted ORFs include three PPE family genes (PPE16, PPE38, and PPE59), and two genes that encode the cell wall proteins Mmpl12 and LprP. PPE family proteins have been well demonstrated to play a role in the antigenic variation of M. tb strains in humans (Dheenadhayalan et al., 2006), and, therefore, the disruption of these genes might affect the membrane structure and alter the antigen profile of strain 1458 in cattle, thereby probably allowing M. tb to evade the host immune system, facilitating an ongoing infection, and further bacterial transmission. Moreover, seven of these 11 intragenic insertions generated 5-, 4-, 4-, 6-, 3-, 3-, and 4-bp direct repeat sequences next to the IS6110 insertion sites (Table 3).

In accordance to previous studies that determined that intergenic regions are preferential targets for IS6110 insertions, we found that the other 10 IS6110s are intergenic insertions at 10 different sites (Table 4), and that they are located upstream of genes such as phoP, esxS, dnaN, ctpD, the PE22-encoding gene, four hypothetical protein-encoding genes, and one alanine-rich protein-encoding gene. These genes encode proteins that play significant roles in M. tb infection and the development of the host immune response. For example, phoP encodes one component of the PhoP-PhoR two-component signaling system, which orchestrates a complex transcriptional program that is essential for mycobacterial growth and virulence (Soto et al., 2004); esxS encodes a member of the 6-kDa early secretory antigenic target (ESAT-6) family of proteins, which are immunogenic and protective (Meher et al., 2006); dnaN encodes a key component of DNA polymerase III, which is a protein essential for many important DNA transactions including replication and repair (Machaba et al., 2016); ctpD encodes a member of the Co²⁺/Ni²⁺-transporting P_1B4-ATPase sub-group, which is important for Co²⁺ and Ni²⁺ homeostasis (Raimunda et al., 2014) and PE22 belongs to the PE/PPE family, which plays an important role in antigenic variation, as mentioned above (Dheenadhayalan et al., 2006). In this study, the 10 intergenic IS6110 insertions are located in the upstream regions of genes, four of which are in the same orientation (5′ to 3′) as the associated genes, while the other six had the reverse orientation. In addition, six of the 10 intergenic insertions generated direct repeat sequences of 3-, 4-, 3-, 4-, 5-, and 3-bp (Table 4).

It was previously demonstrated that IS6110 can act as a promoter when it is inserted in the upstream region of a gene near the transcriptional start site. Firstly, with PCR, we verified IS6110 insertions in the regions of four typical genes. As expected, the amplification products including IS6110 upstream insertion and the genes esxS, phoP, dnaN, and ctpD were 1,784, 1,740, 1,739, and 1,702 bp, respectively. In contrast, for both M. tb H37Rv and M. bovis AF2122_97, because IS6110s were absent in the upstream of these four genes, PCR products are in size of 429, 387, 386, and 349 bp correspondingly (Figure 4).

Then, the copy numbers of esxS, phoP, dnaN, and ctpD transcripts were determined by absolute qPCR. First of all, the transcription of all four genes in M. tb 1458 was significantly higher than M. tb H37Rv (p < 0.01), although the absolute copies of transcripts of the four genes are different; and when the four gene transcription of these three strains was compared, the M. bovis strain expressed significantly highest activity, M. tb 1458 followed, and M. tb H37Rv demonstrated the least activity (Figure 5).

A search for orthologous genes in the genomes of M. tb strain 1458 (GenBank accession no. CP013475), M. tb strain H37Rv (GenBank accession no. NC_000962), and M. bovis strain AF2122_97 (GenBank accession no. NC_002945) identified various clusters of orthologs that were specific to each genome or shared by two or three of the genomes (Figure 6, Table S3).

A total of 3,562 orthologs were identified as the core genome of the three genomes (Figure 6, Table S3-1). These orthologous genes account for 82.8% of the total genes of the three genomes (identified as the pan genome). Among them, 2,818 orthologous genes belong to the COGs classes that are essential for mycobacterial survival and nutrition utilization, including the categories of translation, ribosomal structure, and biogenesis, DNA replication, recombination and repair systems, and carbohydrate, amino acid, nucleotide, and inorganic ion transport and metabolism systems, as well as transcription and energy production and conversion systems. In contrast, 744 of the orthologous genes did not belong to any COGs classes, and they included genes encoding 322 hypothetical proteins, 172 virulence-related factors, 58 transmembrane proteins, 51 membrane proteins, 21 integral membrane proteins, and 10 transcriptional regulatory proteins. The virulence-related factors encoded by the core genome include 127 PE/PE-polymorphic GC-rich repetitive sequence (PGRS)/PPE family proteins, 65 lipoproteins, 47 antitoxins, 42 exported/secreted proteins, 32 epithelial-restricted with serine box (ESX) secretion system proteins (Meher et al., 2006), 23 mammalian cell entry (MCE)-associated proteins (Zhang et al., 2012), 13 ESAT-6 like proteins, 11 toxins (Ramage et al., 2009), four low-molecular-weight antigens, two immunogenic proteins, one integration host factor, and other related factors.

Unique genes were also identified in the three genomes (Figure 6, Tables S3-2–S3-4). One hundred and forty-two genes were considered to be specific for M. tb strain 1458 (Figure 6, Table S3-2). Among them, 70% (100/142) genes were identified as encoding hypothetical proteins, including 37 truncated genes and 63 hypothetical protein-encoding that are unique to the strain 1458 genome, compared with the other two genomes. In addition to these hypothetical protein-encoding, 26 other functional genes were also identified as being specific to the strain 1458 genome, including 13 PE-PGRS/PPE family-encoding genes, eight transposase-encoding genes, and one putative ESAT-6 like protein-encoding gene, esxL (BTB1458_4012). The remaining 16 specific genes are truncated genes, accounting for 37.3% of the specific genes. Most of the truncated genes have orthologs in the other two genomes (Figure 6, Tables S3-5–S3-7). Thirty-nine genes were considered to be specific for M. tb strain H37Rv (Figure 6, Table S3-3), including 14 genes that encode hypothetical proteins, 11 PE-PGRS/PPE family protein-encoding genes, four antitoxin-encoding genes, and two ESAT-6-like protein-encoding genes. Ninety-eight genes were considered to be specific for M. bovis strain AF2122_97 (Figure 6, Table S3-4), including 51 genes that encode hypothetical proteins, six PE/PE-PGRS/PPE family protein-encoding genes, four transmembrane transport protein-encoding genes, four transposase-encoding genes, and one lipoprotein-encoding gene, lppOa.

Overall, the specific hypothetical proteins and PE-PGRS/PPE family proteins seem to be more responsible for the different virulence characteristics of the three strains. Because half of the genes encoding the hypothetical proteins and most of the genes encoding the PE-PGRS/PPE family proteins are orthologous, the differences at the protein level may be attributed to frameshift mutations in these genes. This coincides with the highly similar structures of these three genomes (Figure S1).

A search for SNPs between the M. tb 1458 and H37Rv genomes identified 2,439 SNPs, including 1,400 non-synonymous SNPs, 838 synonymous SNPs, and 201 intergenic SNPs (Table S4-1). The non-synonymous SNPs affected 715 genes in M. tb strain H37Rv. However, the comparison between the M. tb 1458 and M. bovis AF2122_97 genomes identified more SNPs (2,958), including 1636 non-synonymous SNPs, 980 synonymous SNPs, and 342 intergenic SNPs. The non-synonymous SNPs affected 1,037 genes in M. bovis strain AF2122_97 (Table S4-2).

Because the 181 non-synonymous SNP-affected genes in M. tb strain H37Rv and the 246 non-synonymous SNP-affected genes in M. bovis strain AF2122_97 were involved in different pathways, a KEGG pathway enrichment analysis was performed, and it revealed significant enrichment of these genes in some pathways (Tables 5, 6). In the M. tb H37Rv genome, the non-synonymous SNP-affected genes were enriched in carbohydrate metabolism, biosynthesis of other secondary metabolites, and replication and repair (Table 5), while in the M. bovis AF2122_97 genome, the genes were enriched in carbohydrate metabolism and transcription (Table 6). Among them, there is one non-synonymous SNP-affected gene (aceAa) (Rv1915) that encodes an isocitrate lyase that is specific for the M. tb H37Rv genome. It is known that isocitrate lyase mediates the broad antibiotic tolerance of M. tb (Nandakumar et al., 2014), and two isocitrate lyase genes, namely aceAa and aceAb, were found to commonly exist in these three genomes. However, the aceAa gene, which has a single-base-pair overlap resulting in two ORFs (aceAa/b), is specific for the M. tb H37Rv genome. However, it is read as one continuous ORF in M. bovis strain AF2122_97 and M. tb strain 1458, as well as other mycobacterial strains (Zu Bentrup et al., 1999; Nandakumar et al., 2014; Shukla et al., 2015). We speculate that the disruption of aceA in M. tb strain H37Rv may be the result of long-term laboratory growth in which lipids were not the main carbon source, while the redundancy in isocitrate lyase may provide a greater chance for bacterial survival in host cell debris where lipids are used as a carbon source, indicating that the aceA gene may be important in host selection (Zu Bentrup et al., 1999; Shukla et al., 2015).

A previous study determined that M. tb strain 1458, which is of bovine origin, belongs to the M. tb Beijing family genotype (Chen et al., 2009). Therefore, we compared the whole genome sequence of M. tb strain 1458 with those of three other Beijing family strains that were isolated from human TB patients, including two Chinese Beijing genotype strains [M. tb CCDC5079 (NC_021251) and CCDC5180 (NC_017522)] and one Indian strain [Beijing/NITR203 (NC_021054)]. Basically, M. tb strain 1458 is most closely related to the CCDC5079 strain, followed by the CCDC5180 strain (Figure 3). Coincidently, this variation of genetic distance agreed with the drug susceptibility phenotype, as the CCDC5079 strain is sensitive to all four first-line drugs, while the CCDC5180 strain is resistant to all four drugs (Zhang et al., 2011).

A pan-genome analysis and a Mauve alignment were used to compare the strains. The pan-genome analysis identified 3,373 orthologs as the core genome of the four genomes. These orthologs account for 74.6% of the 4,523 total genes of the four genomes (identified as the pan genome) (Table S5-1). In the group of 502 orthologs that are shared by any of the three genomes, 93 orthologs are shared by strains 1458, CCDC5180, and CCDC5079, and these include 15 PE/PPE family proteins and two multidrug resistance efflux proteins (BTB1458_1408 and BTB1458_1409) that were only absent in the M. tb Beijing/NITR203 genome (Table S5-2). Meanwhile, 395 orthologs are shared by the strain 1458, CCDC5079, and NITR203 genomes, and these include 50 antitoxins that are common to the three genomes, while only two antitoxins were found in the core genome (Tables S5-1, S5-2). The toxin-antitoxin (TA) system of these four genomes mainly consists of VapBC family proteins. VapC toxins, specifically their PilT N-terminus (PIN) domains, act as ribonucleases that cleave RNA molecules, thereby reducing the rate of translation (Arcus et al., 2011). Fifty-one toxins, including 28 VapC toxins, 15 putative toxins, two MazF toxins, and one ParE2 toxin were found in the core genome. It has been reported that TA systems are important for M. tb pathogenesis (Ramage et al., 2009); thus, the lack of antitoxins may indicate some loss of function in the CCDC5079 and CCDC5180 genomes. A total of 201 orthologs that are shared by any two genomes were identified (Table S5-3). Most of these orthologs are hypothetical proteins of unknown function.

Unique genes were also identified in each genome (Table S5-4). For the strain 1458, 33 genes are considered to be unique, including 30 hypothetical protein-encoding genes, one PE-PGRS family protein-encoding gene, and one putative collagen-like protein-encoding gene. For M. tb strain CCDC5079, 70 genes are considered to be unique, including 57 hypothetical protein-encoding genes, two lipoprotein-encoding genes, one PE-PGRS family protein-encoding gene, one transposase-encoding gene, and one putative collagen-like protein-encoding gene. For the CCDC5180 genome, 43 genes are considered to be unique, including 22 hypothetical protein-encoding genes, 10 PE/PPE family protein-encoding genes, and two MCE family protein-encoding genes. For the NITR203 genome, 183 genes are considered to be unique, including 110 hypothetical protein-encoding genes, 18 PE/PPE family protein-encoding genes, nine transposase-encoding genes, eight phiRv1 phage-related protein-encoding genes, one phospholipase C4 (MTP40 antigen)-encoding gene, one putative membrane-associated phospholipase C2-encoding gene, one polyketide synthase (Pks13)-encoding gene, one putative ESAT-6-like protein-encoding gene, one putative toxin (VapC16)-encoding gene, and one antitoxin-encoding gene.

In summary, the unique genes of the four genomes totaled 329, and among them, the most common category is that of hypothetical genes, comprising 66.5% of the unique genes, followed by PE/PPE family protein-encoding genes, which comprise 9.1% of the unique genes (Table S5-4). For the individual genomes, strain 1458 has the fewest unique genes (33), 91% (30) of which are hypothetical genes, and 3% (1) of which are PE/PPE family protein-encoding genes, while strain NITR203 has the most unique genes (183), mainly resulting from three inserted fragments (described below), 60.1% (110) of which are hypothetical genes, and 9.8% (18) of which are PE/PPE family protein-encoding genes.

Through the Mauve alignment, three unique insertion fragments were identified in the NITR203 genome (Figure S2). The first (nt 79575 to 83038) affects four genes, including those encoding RNA-directed DNA polymerase, a hypothetical protein, a glutamine-transport transmembrane protein ATP-binding cassette transporter, and a glutamine-transport ATP-binding cassette transporter (nt 479312910 to 479312913). The second (nt 1779126 to 1788367) affects 13 genes that consist of a phiRv1 phage-related gene cluster (nt 479314486 to 479314498). The third (nt 1988260 to 1998428) affects nine genes, including those encoding three transposases, one cutinase, Cut1, one PE-PGRS family protein, one acyltransferase, and three hypothetical proteins (nt 479314672 to 479314680). These genes were also identified as unique in the pan-genome analysis (Table S5-4).

Comparison of the M. tb 1458 genome with three Beijing genomes (CCDC5079, CCDC5180, and NITR203) of human origin revealed that 315 SNPs are specific to CCDC5079, 1,005 SNPs are specific to CCDC5180, and 2,762 SNPs are specific to NITR, while 133 SNPs are unique to the strain 1458 (Table S6). These unique SNPs of the strain 1458 caused non-synonymous variations in 70 genes. A KEGG pathway enrichment analysis only revealed enrichment in RNA polymerase (p < 0.01). Moreover, 127 of the 133 unique SNPs in strain 1458 are either different to those in the M. bovis genome, implying that strain 1458 might have evolved a virulence mechanism that differs from those of M. bovis or the other M. tb Beijing family strains of human origin.

As described above, M. tb strain 1458 has 21 copies of IS6110, 11 of which are distributed within ORFs, while 10 lie between ORFs. Thus, we analyzed the distribution of IS6110 in two other Chinese Beijing genotype strains of human origin, excluding the Indian strain NITR203. The results showed that 18 and 21 copies of IS6110 are located in the CCDC5180 and CCDC5079 genomes, respectively (Table S7).

Compared with M. tb strain 1458, the strain CCDC5079 has two special insertions in a cation efflux system protein-encoding gene and hsdM2, in addition to the copies that are present in hypothetical protein-, PE/PPE family protein-, and transposase-encoding genes. However, neither strain was shown to contain an insertion at the upstream of phoP and esxS genes, suggesting that the virulence of strain 1458 may differ from that of the other two strains.

We analyzed the indels between M. tb strain 1458 and two other Chinese Beijing genotype strains, and the unique insertions of each strain are shown in Table S8. There are 54 unique insertion fragments in these three strains, and they are mainly enriched in genes encoding PE/PPE family proteins (16/54, 29.6%), transposases (17/54, 31.5%), and hypothetical proteins (9/54, 16.6%). For the strain 1458, 15 unique insertion fragments were identified. Of these, a 292 bp insertion in the 5′ terminus of vapB16 is considered to be unique. The vapB16 is a member of the TA loci, and this locus is conserved in other pathogenic strains of the MTBC, including M. bovis and M. avium. For strain CCDC5180, 22 unique insertion fragments were identified. Of these, two insertion fragments, 288 and 302 bp, respectively, were identified in protoheme IX farnesyltransferase. For strain CCDC5079, 17 unique insertion fragments were identified, and of these, a 1,731 bp insertion fragment, involving 542 bp at the 3′ terminus of a gene encoding a member of the phthalate permease family, was unique.

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.