The ability to genetically differentiate among microbial strains facilitates tracing the origins and spread of bacterial pathogens, including Mycobacterium tuberculosis (Mtb), the causative agent of tuberculosis (TB). Various methods for genetically typing clinical strains of Mtb have been developed [1-4]. This includes typing strains through copy number variation (CNV) in mycobacterial interspersed repetitive units (MIRUs), which are minisatellite loci in the Mtb genome [5]. This approach, termed MIRU-VNTR (variable number of tandem repeats; hereafter, MIRU) typing, distinguishes genetically divergent strains rapidly and with relatively high accuracy [6]. As a result, MIRU analysis has been employed in a wide array of epidemiological studies [7-10], where an identical MIRU profile between isolates is interpreted as a recent transmission event [11]. This is based on the assumption that over short periods of time, a change in MIRU copy number is unlikely. Therefore, accurate estimation of the MIRU mutation rate is essential to infer the relationship between transmitted strains based on typing profile. The limits of this approach are highlighted by a recent study of a TB outbreak in British Columbia, where MIRU typing identified a single clonal outbreak, yet higher resolution WGS established two separate, simultaneous outbreaks [12].

Approximating the rate of MIRU CNV based on the profiles of clinical strains has proven to be a challenge. The approaches taken by various groups to estimate a rate have relied primarily on modeling of epidemiological data and meta-analyses [13-16]. The discrepancies in methodology, loci analyzed, and underlying assumptions between different studies have resulted in estimates ranging over two orders of magnitude. For example, Grant et al. (2008) analyzed copy number changes between Mtb lineages and used previous research estimating the time of the most recent common ancestor between lineages to estimate a MIRU mutation rate, yielding a per locus, per year rate of 1.05 × 10^-5. Reyes and Tanaka (2010) used an infinite alleles model to define a relative rate and then benchmarked this rate against estimates of the IS6110 mutation rate to infer a MIRU mutation rate (7.00 × 10^-4- 1.50 × 10^-2). More recently, Aandahl et al. (2012) developed a stepwise mutation model for MIRU evolution and utilized Bayesian statistics to estimate the MIRU mutation rate of previously gathered epidemiological data (3.55 × 10^-3), supporting the estimates provided by Reyes and Tanaka (2010).

Here, we seek to experimentally determine a mutation rate by assessing CNV in Mtb strains isolated from cynomolgus macaques [17], an animal model of Mtb infection that closely recapitulates the course of human disease [18]. The genomes from Mtb isolated from infected macaques were previously sequenced using the Illumina platform [17]. While protocols on detecting single nucleotide polymorphisms (SNPs) and small insertion/deletions (indels) using Illumina sequencing data are well established, it is less clear how to best determine the copy number in minisatellite loci. As WGS becomes increasingly common in epidemiological studies, the ability to establish MIRU copy number from sequencing data becomes important for the analysis of new WGS data in the context of previously existing typing data.

The length of Illumina reads fails to span the majority of MIRU repeats in a locus, which range from one to five or more repeats of 40–100 basepairs each in Mtb. Thus, short read sequencing cannot capture unique sequence and define copy number. Recently, researchers have utilized read depth (RD), a measure of the density of sequencing reads at each nucleotide in the genome, in order to identify CNV [19,20]. This approach relies on the observation that the absolute number of sequencing reads mapped to a reference genome is proportional to the copy number of a particular strain [19]. However, this approach has been primarily utilized to identify large (> 1 kb) variants in human tandem repeats, and RD has not been successfully applied to assess smaller minisatellite variation in microbes. Similarly, attempts to identify structural variation, including CNV, using mate pair distance have been previously employed in human genome studies [21-23]. This approach, termed paired-end mapping (PEM), utilizes the likely distance between paired reads from Illumina sequencing to identify structural variants. When mapping back to a reference genome, if the distance between paired reads is discordant from the expected value, it is suggestive of CNV relative to the reference genome. PEM has been successfully employed to identify large insertions or deletions. However, bacterial minisatellite CNV produces only small discordances relative to the reference genome which fall within the expected distribution of mate pair distance, making such events difficult to detect.

In this study, we sought to define a WGS methodology useful for identifying MIRU CNV by mapping sequencing reads to a single copy-number reference genome. We validated the WGS analyses by resequencing all MIRU loci via Sanger sequencing, which allowed us to assess the accuracy of using RD and PEM approaches to identify minisatellite variations. We then used our WGS and Sanger resequencing data to estimate a MIRU mutation rate during the course of infection. This rate will help guide the analysis of epidemiological data and provide a preliminary understanding of site-specific mutability in Mtb.

Here we have utilized both conventional and next-generation sequencing approaches to define the number of copy variants at MIRU loci that arose during the course of macaque infection. Approaches to detect MIRU CNV by Illumina sequencing are increasingly important given the expanding use of WGS in molecular epidemiology. However, mapping ambiguities due to the repetitive nature of MIRU loci complicate copy number assessment by WGS. In this work, we sought to reduce mapping ambiguity and identify CNV in Mtb by mapping to a single copy-number genome. Sequencing technologies have advanced since the sequencing of these strains, and longer reads with less coverage variability may improve specificity in future work [31,32]. Additionally, it is thought that variability in both RD and PEM arises from local disparities in GC content [25,33]. Future work may exploit advances in correcting for GC bias [19,32,34-36], in combination with the approaches described here, to more reliably detect MIRU CNV from WGS data.

With the MIRU CNV identified in strains isolated from cynomolgus macaques, we have estimated the mutation rate at MIRU loci. Our in vivo MIRU mutation rate, 2.70 × 10^-03 per locus, per year, is most consistent with the highest published estimates. Variability in these estimates may be partially driven by differences in MIRU loci analyzed, the epidemiology of the strains used, and differences in the assumptions of the models used to estimate a rate. The resulting differences in rate estimates have motivated extensive debate in the literature [15,37,38]. Our estimate is derived from an alternative, experimental approach, relying on the cynomolgus macaque model of TB infection to assess the mutation rate of the 24 loci standard used in MIRU typing. Though our analysis is somewhat limited by a relatively small signal, strikingly, our estimate closely aligns with previous rates.

What are the biologic consequences of a high mutation rate at MIRU loci? Interestingly, most MIRUs are located in intergenic regions and are hypothesized to be transcribed as part of a polycistronic operon. Several MIRU elements are located within the coding region of well-described two-component regulatory systems as well as genes essential for virulence and host adaptation [5]. It is interesting to hypothesize that rapid genetic variation at these MIRU loci may have effects on the transcription of the regulon, thus generating population diversity. Some evidence exists that CNV variation at MIRU loci may result in transcriptional changes of the downstream gene within a MIRU locus [39,40], though further characterization is required to establish the biologic relevance of these loci and the role of variation in MIRU elements.

Repeat variation is a well-established means of generating locus specific mutation in other microbial genomes [41,42]. The rate established here is comparable to tandem repeat variation rates in other organisms. For example, in Bacillus anthracis the estimated mutation rate of tandem repeat loci is reported to be roughly 10^-05 to 10^-04 per generation [43], while in pathogenic E. coli O157:H7, the rate is reported to be roughly 6.4 × 10^-04 per generation [44]. The MIRU per generation mutation rate is between 5.6 × 10^-06 and 7.64 × 10^-05, using the lower (18 hours) and upper (240 hours) estimates of generation times previously described [17]. It is important to note the loci analyzed in this study were selected from the larger set of MIRU loci for their relative stability, making them ideal for typing. The rate of variation in the remaining loci may vary from the rate reported here, especially in loci previously identified as hypermutable [6]. While further work is needed to determine the biological consequences of MIRU variation, we have shown that there is potential for detecting MIRU variants by WGS and that the rate uncovered from macaque infection is consistent with the highest previous estimates.

CNV in typing markers is an essential tool to differentiate and classify clinical strains, and quantitation of marker variation allows for enhanced interpretation of epidemiological data. In this study, we have used the macaque model of Mtb infection to estimate the MIRU mutation rate during the course of disease, and we have explored the use of WGS to assess MIRU copy number. Subsequent Sanger resequencing confirmed two of the four MIRU indels identified by PEM to a reduced copy genome, and from this we have estimated a per locus, per year mutation rate of 2.70 × 10^-03. This value agrees with the higher published estimations of MIRU mutation rates. Further assessment of RD and PEM as indicators for copy number may streamline minisatellite detection via WGS.

MIRU: Mycobacterial interspersed repetitive unit; Mtb: Mycobacterium tuberculosis; CNV: Copy number variation; WGS: Whole-genome sequencing; RD: Read depth; PEM: Paired-end mapping; TB: Tuberculosis; VNTR: Variable number of tandem repeats; SNP: Single nucleotide polymorphism; Indel: Insertion/deletion; SD: Standard deviation.

The authors declare that they have no competing interests.

MNR and CBF designed and performed molecular studies, conducted the data analyses, prepared the figures and drafted the manuscript. MRC analyzed sequence data. PLL and JLF conducted the infection of the cynomolgus macaques, determined clinical state and acquired bacterial strains on necropsy. SMF designed the study, supervised experimental and molecular studies, and drafted the manuscript. All authors edited the manuscript.