This study assembles DNA adenine methylomes for 93 Mycobacterium tuberculosis complex (MTBC) isolates from seven lineages paired with fully-annotated, finished, de novo assembled genomes. Integrative analysis yielded four key results. First, methyltransferase allele-methylome mapping corrected methyltransferase variant effects previously obscured by reference-based variant calling. Second, heterogeneity analysis of partially active methyltransferase alleles revealed that intracellular stochastic methylation generates a mosaic of methylomes within isogenic cultures, which we formalize as ‘intercellular mosaic methylation’ (IMM). Mutation-driven IMM was nearly ubiquitous in the globally prominent Beijing sublineage. Third, promoter methylation is widespread and associated with differential expression in the ΔhsdM transcriptome, suggesting promoter HsdM-methylation directly influences transcription. Finally, comparative and functional analyses identified 351 sites hypervariable across isolates and numerous putative regulatory interactions. This multi-omic integration revealed features of methylomic variability in clinical isolates and provides a rational basis for hypothesizing the functions of DNA adenine methylation in MTBC physiology and adaptive evolution.

Whole-genome sequencing (WGS) is fundamental to Mycobacterium tuberculosis basic research and many clinical applications. Coverage across Illumina-sequenced M. tuberculosis genomes is known to vary with sequence context, but this bias is poorly characterized. Here, through a novel application of phylogenomics that distinguishes genuine coverage bias from deletions, we discern Illumina ‘blind spots’ in the M. tuberculosis reference genome for seven sequencing workflows. We find blind spots to be widespread, affecting 529 genes, and provide their exact coordinates, enabling salvage of unaffected regions. Fifty-seven pe/ppe genes (the primary families assumed to exhibit Illumina bias) lack blind spots entirely, while the remaining pe/ppe genes account for 55.1 % of blind spots. Surprisingly, we find coverage bias persists in homopolymers as short as 6 bp, shorter tracts than previously reported. While G+C-rich regions challenge all Illumina sequencing workflows, a modified Nextera library preparation that amplifies DNA with a high-fidelity polymerase markedly attenuates coverage bias in G+C-rich and homopolymeric sequences, expanding the ‘Illumina-sequenceable’ genome. Through these findings, and by defining workflow-specific exclusion criteria, we spotlight effective strategies for handling bias in M. tuberculosis Illumina WGS. This empirical analysis framework may be used to systematically evaluate coverage bias in other species using existing sequencing data.

BackgroundThe genetic basis of virulence in Mycobacterium tuberculosis has been investigated through genome comparisons of virulent (H37Rv) and attenuated (H37Ra) sister strains. Such analysis, however, relies heavily on the accuracy of the sequences. While the H37Rv reference genome has had several corrections to date, that of H37Ra is unmodified since its original publication.ResultsHere, we report the assembly and finishing of the H37Ra genome from single-molecule, real-time (SMRT) sequencing. Our assembly reveals that the number of H37Ra-specific variants is less than half of what the Sanger-based H37Ra reference sequence indicates, undermining and, in some cases, invalidating the conclusions of several studies. PE_PPE family genes, which are intractable to commonly-used sequencing platforms because of their repetitive and GC-rich nature, are overrepresented in the set of genes in which all reported H37Ra-specific variants are contradicted. Further, one of the sequencing errors in H37Ra masks a true variant in common with the clinical strain CDC1551 which, when considered in the context of previous work, corresponds to a sequencing error in the H37Rv reference genome.ConclusionsOur results constrain the set of genomic differences possibly affecting virulence by more than half, which focuses laboratory investigation on pertinent targets and demonstrates the power of SMRT sequencing for producing high-quality reference genomes.Electronic supplementary materialThe online version of this article (doi:10.1186/s12864-017-3687-5) contains supplementary material, which is available to authorized users.

For decades, pyrazinamide has served as a cornerstone of tuberculosis therapy. Unlike any other antitubercular drug, pyrazinamide requires an acidic environment to exert its action.

Pyrazinamide (PZA) is a widely used antitubercular chemotherapeutic. Typically, PZA resistance (PZA-R) emerges in M. tuberculosis strains with existing resistance to isoniazid and rifampicin (MDR) and is conferred by loss-of-function pncA mutations that inhibit conversion to its active form, Pyrazinoic acid (POA). PZA-R departing from this canonical scenario is poorly understood. Here, we genotype pncA and purported alternative PZA-R genes (panD, rpsA, and clpC1) with long-read sequencing of nineteen phenotypically PZA mono-resistant isolates collected in Sweden and compare their phylogenetic and genomic characteristics to a large set of MDR PZA-R (MDRPZA-R) isolates. We report the first association of ClpC1 mutations with PZA-R in clinical isolates, in the ClpC1 promoter (clpC1p-138) and N-terminal (ClpC1Val63Ala). Mutations have emerged in both these regions under POA selection in vitro and ClpC1N-terminal has been implicated further, through its POA-dependent efficacy in PanD proteolysis. ClpC1Val63Ala mutants spanned 4 Indo-oceanic sublineages. Indo-oceanic isolates invariably harbored ClpC1Val63Ala and were starkly overrepresented (OR=22.2, p <0.00001) among PZA mono-resistant isolates (11/19) compared to MDRPZA-R isolates (5/80). The genetic basis of Indo-oceanic isolates’ overrepresentation in PZA mono-resistant TB remains undetermined, but substantial circumstantial evidence suggests ClpC1Val63Ala confers low-level PZA resistance. Our findings highlight ClpC1 as potentially clinically relevant for PZA-R and reinforce the importance of genetic background in the trajectory of resistance development.

Mycobacterium tuberculosis , the primary causative agent of tuberculosis, kills more humans than any other infectious bacterium. Yet 40% of its genome is functionally uncharacterized, leaving much about the genetic basis of its resistance to antibiotics, capacity to withstand host immunity, and basic metabolism yet undiscovered.

The genetic basis of virulence in Mycobacterium tuberculosis has been investigated through genome comparisons of its virulent (H37Rv) and attenuated (H37Ra) sister strains. Such analysis, however, relies heavily on the accuracy of the sequences. While the H37Rv reference genome has had several corrections to date, that of H37Ra is unmodified since its original publication. Here, we report the assembly and finishing of the H37Ra genome from single-molecule, real-time (SMRT) sequencing. Our assembly reveals that the number of H37Ra-specific variants is less than half of what the Sanger-based H37Ra reference sequence indicates, undermining and, in some cases, invalidating the conclusions of several studies. PE PPE family genes, which are intractable to commonly-used sequencing platforms because of their repetitive and GC-rich nature, are overrepresented in the set of genes in which all reported H37Ra-specific variants are contradicted. We discuss how our results change the picture of virulence attenuation and the power of SMRT sequencing for producing high-quality reference genomes. * Equal contribution † Corresponding author: faramarz@sdsu.edu 1 . CC-BY-ND 4.0 International license peer-reviewed) is the author/funder. It is made available under aThe copyright holder for this preprint (which was not . http://dx.doi.org/10.1101/064840 doi: bioRxiv preprint first posted online Jul. 19, 2016; Tuberculosis is a serious and pervasive public health problem [1]. It is a disease 1 caused by infection of bacteria from the Mycobacterium tuberculosis complex 2 (MTBC). The reference strain, Mycobacterium tuberculosis H37Rv, has an at-3 tenuated counterpart known as H37Ra that is available for studies where facil-4 ities to handle virulent samples are lacking. H37Ra exhibits a distinct colony 5 morphology, an absence of cord formation, decreased resistance to stress and 6 hypoxia, and attenuated virulence in mammalian models [2][3][4]. The H37Ra 22In this study, we sequenced and assembled the genome of M. [14, 15]. This insertion was the heterogeneous inser-56 tion responsible for the discrepant contig ends in our raw genome assembly. 57Such heterogeneity implies either a lack of selection pressure on the insertion in 58 culture, a recent emergence of the insertion, or both. 59The 3456bp insertion in ppe54 with respect to H37RaJH incidentally corre-60 sponds to a tandem duplication of a 1728bp sequence at the same site in H37Rv H37Rv and CDC1551 as "H37Ra-specific". These mutations fall within or ad-95 jacent to (which we term "affecting") 56 genes in H37Rv, which we refer to as The copyright holder for this preprint (which was not . http://dx.doi.org/10.1101/064840 doi: bioRxiv preprint first posted online Jul. 19, 2016; colleagues also discovered sequencing errors in the H37Rv reference sequence [5], 98 a number of which were corrected in NC 000962.3 [9], the version used in our 99 study. 100To see how well the HC genes are supported by our assembly of H37Ra, we 101 determined variants with respect to H37Rv ...

Whole genome sequencing (WGS) is fundamental to M. tuberculosis basic research and many clinical applications. Coverage across Illumina-sequenced M. tuberculosis genomes is known to vary with sequence context, but this bias is poorly characterized. Here, through a novel application of phylogenomics that distinguishes genuine coverage bias from deletions, we discern Illumina "blind spots" in the M. tuberculosis reference genome for seven sequencing workflows. We find blind spots to be widespread, affecting 529 genes, and provide their exact coordinates, enabling salvage of unaffected regions. Fifty-seven PE/PPE genes (the primary families assumed to exhibit Illumina bias) lack blind spots entirely, while remaining PE/PPE genes account for 55.1% of blind spots. Surprisingly, we find coverage bias persists in homopolymers as short as 6 bp, shorter tracts than previously reported. While GC-rich regions challenge all Illumina sequencing workflows, a modified Nextera library preparation that amplifies DNA with a high-fidelity polymerase markedly attenuates coverage bias in GC-rich and homopolymeric sequences, expanding the "Illumina-sequencable" genome. Through these findings, and by defining workflow-specific exclusion criteria, we spotlight effective strategies for handling bias in M. tuberculosis Illumina WGS. This empirical analysis framework may be used to systematically evaluate coverage bias in other species using existing sequencing data.