BACKGROUND
in spite of effective chemotherapy against tuberculosis, this disease is still a global health problem and a leading cause of death worldwide  <cit> .
the causative organism, mycobacterium tuberculosis, is an intracellular pathogen that has infected humans since ancient times  <cit> .
tuberculosis disease  results from an intricate interaction between the host immune system’s efforts to control the infection and the pathogen’s ability to grow and persist within the host.
thus, infection with the tubercle bacillus has variable outcomes that range from sterilizing immunity to active tb  <cit> .
active disease occurs in only 5– <dig> percent of immunocompetent individuals, with pulmonary tuberculosis being responsible for transmission in the community  <cit> .
in most cases, however, the infection is controlled by the host’s immune system and can lead to the establishment of a long-term latent infection which can reactivate later in life  <cit> .
there is no evidence of recent genetic exchange in m. tuberculosis which is thus considered to have a clonal population structure, with strains having almost identical nucleotide sequences  <cit> .
however, there is substantial intra-species genetic diversity that can affect disease outcome  <cit> .
the insertion element is <dig> is an important source of genomic variability for m. tuberculosis and is distributed in different positions across the bacterial chromosome  <cit> .
the copy number of this element is highly variable, but most strains contain between  <dig> and  <dig> instances.
as a result, is <dig> has been widely used as an epidemiological marker in tuberculosis by dna fingerprinting using restriction fragment length polymorphism   <cit> .
these insertions also represent “in vivo” transposition events that provide information regarding genes required for human infection and disease.
in addition, the persistence of this element in the m. tuberculosis genome raises the possibility that it might also drive phenotypic variability or affect strain fitness.
it has been demonstrated that is  <dig> transposition can be stimulated by specific stress conditions and thus contribute to genetic diversity in circulating m. tuberculosis clinical isolates  <cit> .
is <dig> insertions can affect gene expression by interrupting protein-coding genes, by mediating recombination events that result in deletions and inversions, and by up-regulating the expression of nearby genes due to a promoter located within the transposable element  <cit> .
thus, although is <dig> is primarily thought of as a valuable epidemiological tool, its prevalence and effect on genome function prompted us to take a deeper look at the distribution and patterns of is  <dig> insertions by conducting a broad survey of circulating strains.
a previous examination of is <dig> insertion sites in  <dig> m. tuberculosis isolates, which was accomplished by pcr amplification of the target sequences followed by cloning and sequencing, identified  <dig> unique insertions sites in only  <dig> genes, suggesting that a large gene repertoire is required for human infection  <cit> .
however, the technical difficulties and biases due to the amplification and cloning associated with this approach prevented the analysis of a large number of isolates and the detection of all of the insertions present in each strain.
to overcome these limitations and assess the extent of naturally occurring is <dig> insertions during infection and transmission, we used illumina® sequencing technology to identify is  <dig> flanking regions in more than  <dig> m. tuberculosis isolates from representative collections of clinical isolates from europe and south america.
this high-throughput approach obviates the need for cloning and allows analysis of many strains in a cost-effective manner.
we were able to identify approximately  <dig>  is  <dig> insertion sites that interrupt almost  <dig> genes, which together represent <10% of the m. tuberculosis genome.
thus, we substantially increased the survey of naturally occurring in vivo is  <dig> transpositions, providing new insights regarding their distribution, diversity and role in strain classification
RESULTS
insertion site sequencing
we developed a method for the simultaneous identification of insertion sites  in up to  <dig> different strains by barcoding samples and using illumina® sequencing technology to sequence the genomic dna flanking both ends of an insertion sequence .
we selected  <dig> m. tuberculosis strains for analysis of the is  <dig> insertion sites .
this set included the reference strain m. tuberculosis h37rv  and six other strains  that were used to validate the method by verification of insertion sites using pcr .
all of the strains included in this study were previously characterized by is <dig> rflp and spligotyping.
*ppv: positive predicted value.
for each reference strain the number of rflp observed bands is reported as well as the number of is sites identified using is-seq by one or both termini of the insertion and the number of confirmed sites by pcr.
sensitivity and positive predicted value is calculated as indicated in methods.
to generate a sequencing library we randomly sheared genomic dna and ligated barcoded adapters to the fragment ends.
fragments that include the terminal end of the is <dig> were then amplified in a nested pcr using primers that incorporate additional barcodes.
finally, the amplified products were pooled, purified and sequenced.
the resulting fragments contain either the 3’ or 5’ terminal sequence of is  <dig>  the flanking genomic sequence, and two different barcodes at each terminus.
since most of the enzymatic steps are “add only” reactions, there is no material loss due to purification steps and the protocol can be easily automated and performed in 96-well plates.
furthermore, reactions are pooled after both barcodes have been added, so that the gel purification of hundreds of libraries can be performed using a single gel lane.
using this protocol, we obtained more than  <dig> million “paired-end” sequences using an illumina® ga ii machine.
the sequence reads contained an even distribution of barcodes , showing that pooling before the gel extraction without individual quantification introduced no bias on the distribution.
we were able to identify sequences for all of the  <dig> strains with an average of  <dig>  reads per strain.
to increase the precision of mapping, the short reads were sorted by sample and assembled to generate contigs of up to  <dig> nucleotides before being mapped to the reference m. tuberculosis genome.
all assembled contigs were mapped to at least one of the available genomes.
the sequences that did not assemble in contigs were further mapped as single reads to the different genomes and less than 5% of these reads could not be mapped.
based on an estimate of the average sequence coverage obtained per is  <dig> element , a cutoff point was defined  to remove strains with low coverage, which are likely to produce unreliable results .
after this filtering,  <dig> strains were retained for analysis .
we included a negative control species in this study that lacks is <dig>  and obtained only few and scattered reads without significant mapping.
these results demonstrate the efficiency and specificity of the protocol and primer design.
accuracy of the method
the accuracy of is-seq was assessed using seven m. tuberculosis strains for which  <dig> insertions had been previously characterized by rflp fingerprinting.
is-seq identified these  <dig>  as well as  <dig> additional insertions.
of these additional insertions, all but two were confirmed by pcr .
these two false positives corresponded to two insertions in the h37rv reference strain  that were mapped to only one terminus of the is <dig> element.
thus, by requiring sequence from both ends of the is  <dig> element 100% specificity could be achieved, at the cost of sensitivity .
the identification of additional insertion sites using is-seq was not altogether surprising given the limited resolution of the rflp fingerprinting method.
based on the analysis of these strains, the is-seq method had a sensitivity of 100% and a positive predictive value of 98% .
it was also interesting to find that while the col-h37rv strain used here had  <dig> copies of is <dig>  the sequenced genome of m. tuberculosis h37rv  had  <dig> is  <dig> elements, with an extra copy being part of a doublet in positions  <dig> , <dig> and  <dig> , <dig>  <cit> , highlighting differences in the genetic makeup of these reference strains.
the reproducibility of the method was assessed by including two sets of duplicate samples of the same dna ; the agreement on the number and insertion position within each pair of dnas was greater than 90% .
a hallmark of the is <dig> element insertion is the generation of a 3- <dig> bp duplication at the insertion site  <cit> .
we expected to identify this duplication in those cases where we were able to map the insertion site of the is <dig> element on both termini .
in  <dig> % of these cases we identified a  <dig> bp  or  <dig> bp  duplication.
taken together, these results demonstrate that the is-seq method accurately and reproducibly identifies and maps is  <dig> insertion sites.
genomic distribution of the is <dig> element
the is-seq method identified  <dig>  insertion sites in the  <dig> strains analyzed .
these represented  <dig> different insertions that affected  <dig> loci , 66% of which were in coding regions  and 34% were intergenic .
this represents a significant bias towards non-coding regions, since intergenic regions account for only ~10% of the m. tuberculosis genome.
an analysis of the distribution and frequency of independent insertion sites throughout the genome also showed a significant location bias.
a poisson test for the distribution of insertion sites showed that there were forty-four  <dig> kb-windows  with significantly more insertions than expected .
when  <dig> kb windows were used, regions with significantly more or less insertions were also identified .
an analysis of these regions revealed no significant correlation with gc content .
we found multiple insertion sites in 40% of the affected loci .
of these loci,  <dig> contained  <dig> or more independent sites and thus were labeled as “hotspots”, with the “hottest” locus containing  <dig> different insertion elements .
several of these “hotspots” have been previously described, such as mmpl <dig> , which is thought to be involved in fatty acid transport, or plcd, which encodes a phospholipase c, and intergenic regions like dnaa-dnan or the region between genes rv2813- <dig> .
this latter locus corresponds to a cluster of regularly interspaced short palindromic repeats  known as the direct repeat  region, which is the template for the spoligotyping technique  <cit> .
we identified  <dig> insertions in  <dig> different sites within the dr region , one of which was present in  <dig>  of the evaluated strains, including the reference h37rv genome .
interestingly, this is the only is <dig> present in the vast majority of mycobacterium bovis isolates, suggesting that it could have been present before many of the modern lineages diverged from each other.
we therefore refer to this as an ancestral point of insertion.
this insertion was absent in the beijing strains evaluated in our study, which is congruent with the loss of spacers 1– <dig> in the dr region that characterizes this m. tuberculosis family  <cit> .
to test if our analysis was close to saturation in terms of identifying all possible genes that could be naturally interrupted by the is <dig> element, a rarefaction analysis was carried out .
the lack of a plateau indicating saturation in mutable sites shows that additional insertions in coding sequences could still be identified in circulating m. tuberculosis strains.
this study was done with representatives of the major lineages circulating in europe and latin america and falls short of representatives of strains prevalent in africa and asia that may harbor distinct insertion sites.
a regression analysis on the rarefaction curves estimated that between  <dig> and  <dig> genes could support insertion sites, consistent with previous estimates of the number of possible in vivo mutations  <cit> .
the presence of mutable sites results from an ongoing evolutionary process that in m. tuberculosis might not necessarily have reached saturation.
thus, even if we sampled all available isolates there is no guarantee that we could observe such saturation of in vivo insertions.
analysis of disrupted genes
we next characterized the genes disrupted by is <dig> insertions by assigning them to functional categories based on the ncbi clusters of orthologous groups   and the tuberculist database .
disrupted genes were broadly distributed across functional categories, with a significant over-representation of insertions in genes belonging to three categories: 1) pe/ppe, 2) genes in cog from the category “unknown” and 3) cell wall/membrane/envelope biogenesis genes .
the pe and ppe gene families encode glycine-rich proteins that contain repeats of the motifs pro-glu  or pro-pro-glu  and represent about 10% of the predicted proteins in m. tuberculosis <cit> .
within the pe/ppe gene family, the majority of interrupted loci  occurred in ppe genes.
under-represented cog categories included genes involved in lipid and amino acid transport and metabolism, energy production and conversion, as well as the category defined as virulence, detoxification and adaptation in tuberculist .
<dig>  proteins belonging to the pe/ppe families were extracted to a separate category, as they constitute an important family of proteins in m. tuberculosis.
<dig>  the genes for the transposase is <dig> were removed from the category of replication, recombination and repair.
<dig>  cumulative length  of all the genes in a given category.
the probability of under or over-representation of a given functional category is dependent on both the number of genes and their length.
<dig>  represents the number of independent insertion events identified in genes of a given category.
<dig>  probability of over or under-representation of insertion sequences interrupting genes of a given category.
in bold, categories with significant over-, under-representation after bonferroni correction.
bonferroni corrected threshold =  <dig> e- <dig>
<dig>  this category is not part of cog; it is defined in tuberculist .
we next compared the set of orfs harboring an is <dig> element with previous studies that identified m. tuberculosis genes required for growth in vitro or in vivo  <cit> .
twenty-seven genes previously identified as being essential under different experimental conditions were found to contain is <dig> elements .
seven of these genes contained the insertion element in the terminal 20% of the gene, which although less likely to affect function could still do so in proteins with c-terminal active sites or domains.
to rule out the possibility that these observations were due to an artifact of the is-seq method, we validated these insertion sites by pcr.
twenty-four of the  <dig> insertions were confirmed , and one locus failed to amplify, probably due to the limited amount of dna used as starting material .
rv3473c  indicates the position of the insertion within the orf as the percentage of gene upstream of the insertion.
2) insertions verified using pcr: +, pcr positive; -, pcr negative.
is <dig> family signatures
to see how well these is <dig> patterns correlated with known strain classifications, we examined the relationship between is  <dig> distribution and m. tuberculosis strain families .
the use of specific is <dig> insertions as markers for different m. tuberculosis lineages has been described.
for example, the beijing lineages have been described to contain an is  <dig> insertion in the intergenic region between dnaadnan <cit>  and the spoligotype characteristic of the haarlem lineage  is partially due to a characteristic is <dig> in the dr region  <cit> .
despite the limitation of using the is <dig> to accurately discriminate between the different m. tuberculosis phylogenetic lineages  <cit> , we propose that the precise locations of is <dig> insertions could be useful as markers within a lineage and, furthermore, could help to identify previously uncharacterized clades within families.
we therefore used a mutual information metric  to identify the co-occurrence or absence of insertion sites in order to generate a dendrogram using  <dig> of our strains, which belonged to the m. tuberculosis euro-american lineage  and the beijing family , which was also used as an outgroup .
we also included the  <dig> publicly available genomes that belong to either lineage  <dig> or  <dig> .
the tree generated distinguished the main lineage  <dig> families , and subdivisions within families, such as the h <dig>  h <dig> and h <dig> haarlem sub-families defined by spoligotyping.
for the lam family, two main clusters were formed, one composed of lam <dig>  lam <dig>  lam <dig> and the second of lam <dig>  lam <dig>  lam <dig> and tuscany.
the t family also showed sub-clusters that did not reflect the t sub-families described by spoligotyping, but were consistent with findings using snps for classification that suggested that t sub-families are still not well defined  <cit> .
by analyzing the patterns of insertion events in the different families, we identified  <dig> family-specific insertions, consisting of few insertions per family .
these signature is sites helped to understand inconsistencies in spoligotyping .
in particular,  <dig> strains ambiguously classified by spoligotyping as either ‘haarlem orphan, h <dig> and h <dig>  and only as haarlem by rflp, were found to contain a h <dig> subfamily characteristic insertion site and a h <dig> subfamily spoligotype pattern in the dr region, which may have occurred independently as a product of recombination between two is <dig>  interestingly, we were able to identify the insertions in the dr region that are responsible for the spoligotype patterns of the different members of the haarlem family .
two haarlem  <dig> strains were found to be ural strains, consistent with recent reports  <cit> , and the tuscany subfamily was closely related to the lam subfamily rather than the t family , as had been recently suggested  <cit> .
also, the previously described insertion in the dnaa -dnan intergenic region in beijing isolates was identified in all  <dig> strains.
when these family-specific insertions were used to classify the  <dig> strains used in this study, family or subfamily information was provided for  <dig> strains compared to  <dig> that were classified by spoligotyping.
these signatures were also used to classify  <dig> of the  <dig> available complete m. tuberculosis genomes  .
some strains were assigned to specific families, as was the case of the reference h37rv and the h37ra strains that were assigned to the t family, and the f <dig> and haarlem strains that were identified as belonging to lam <dig> and haarlem  <dig>  respectively.
other strains were classified as lam <dig> , lam , x  or beijing ; all in 100% agreement with previously suggested families, when available   <cit> .
the inability to classify the remaining  <dig> cases is due to the fact that they belong to lineages not significantly represented in our study .
these results show that classification of tuberculosis strains by the use of signature is sites is consistent with snp classification and therefore highly reliable.
to address the possibility that insertions could appear independently in different lineages because they confer some selective advantage, such as drug resistance, we performed a conditional mutual information analysis  <cit>  to see if any insertion sites correlated significantly with the epidemiological data collected for each strain .
no association was found between insertion sites present in different phylogenetic families and a specific biological trait, suggesting that it is unlikely that exact, yet independent insertions could occur in different families that confer a fitness advantage .
furthermore, in all cases where a specific insertion sequence was present in multiple strains, the strains containing that insertion belonged to the same phylogenetic lineage, suggesting that occurrence of independent insertions in the same nucleotide position but in different lineages by convergent evolution, although plausible, is very unlikely.
instead, hotspots of insertions occurring in a single locus within few nucleotides of distance were observed .
discussion
the is-seq method was effective at detecting, in a single reaction, multiple is <dig> insertions sites in over  <dig> m. tuberculosis isolates from diverse geographical origins and phylogenetic lineages.
this was achieved using a novel approach with the following features: i) the use of two molecular barcodes per sample to multiplex samples on an illumina® sequencer, allowing the analysis of a combinatorial number of samples  with just  <dig> different barcodes; ii) a high-throughput library construction protocol that uses minimal amounts of starting material, and iii) a novel computational strategy that maximizes the efficiency of sequence identification by assembling reads prior to mapping against the m. tuberculosis sequenced genomes.
this study, combined with the advances on sequencing technologies, provides a highly cost-efficient method for m. tuberculosis strain identification and typing that could be extended to other target amplifications and sequencing.
while the use of  <dig> bp barcodes allows multiplexing of up to  <dig> samples per sequencing reaction,  <dig> bp barcodes may allow pooling up to  <dig>  samples.
additionally, the requirement of just  <dig> reads per insertion site for accurately mapping, combined with new technologies such as miseq that achieve up to  <dig> million reads in just a matter of hours, can revolutionize the field of tuberculosis strain characterization.
is-seq effectively detected loci with multiple independent insertions or insertions in the same location in many strains, and identified new and previously characterized hotspots and insertions that could be considered ancestral.
although an is <dig> element present in the same locus could be due to convergent evolution, which in time could be reflected as insertion hot spots  <cit> , our observations suggest that insertions at the same nucleotide position are more likely due to a common ancestral insertion.
strategies similar to is-seq have been recently developed for identifying microbial transposon insertions sites generated during in vitro mutagenesis  <cit>  and endogenous transposon elements in the human genome  <cit> .
however, our method differs because it can identify multiple insertions sites per strain and simultaneously analyze hundreds of different strains per reaction.
the advantage of is-seq results from a more streamlined library preparation protocol and a more efficient barcoding strategy.
the insertion sites identified in these  <dig> clinical m. tuberculosis isolates significantly extends previous surveys by identifying more that  <dig> new sites and provides, for the first time, a comprehensive overview of the gene repertoire that can be interrupted in vivo by is  <dig>  rarefaction curves indicated that additional genes could still be identified, especially in strains that were not significantly represented in our dataset.
however, the number of loci predicted to sustain in vivo insertions was close to the maximum of about  <dig> genes previously predicted based on an analysis of  <dig> strains  <cit> .
it is also consistent with the observation that the number of naturally occurring in vivo mutations is much lower than the genes predicted to be non-essential using transposon mutagenesis approaches  <cit> , probably due to the detrimental effect of cumulative insertions on in vivo fitness.
the smaller gene set that can be inactivated in vivo when compared to in vitro studies, tied to identification of insertions in previously predicted essential genes, indicates evident differences in the conditions imposed on circulating strains as they adapt to the human host.
in this respect, the overrepresentation of insertions in pe/ppe, cell wall and membrane biogenesis genes and orfs of “unknown function” could reflect adaptation as bacteria alter envelope components in response to pressure imposed by the host environment.
in particular, the presence of is <dig> elements in the pe/ppe gene family reinforces their previously proposed role in immune evasion and is consistent with previous reports  <cit> .
however, the over-representation of certain functional categories could also be due to redundancy in function that might alleviate the effect of gene disruption.
in contrast, underrepresented gene categories, such as energy production and conversion, highlight their importance in processes required for bacterial virulence and survival.
finally, the identification of is <dig> insertions in genes reported to be essential is not implausible, since insertions as well as deletions of some of these genes have already been described  <cit> .
in addition to gene disruption, the is <dig> insertions identified here can also result in transcription up-regulation of downstream genes, as has been reported for an insertion in the promoter region of the two-component system gene phop <cit> .
a similar insertion in the phop promoter was observed here in a fast evolving clade of haarlem  <dig> in latin america .
this survey therefore broadens the list and emphasizes again the difference in requirements between strains grown and studied under laboratory conditions and those subjected to the pressure of in vivo passage and processes of natural infection in the human host.
interestingly, the strain classification based on the distribution of is <dig> using is-seq was congruent with that based on spoligotyping and snps for sequenced genomes.
both is <dig> rflp and spligotyping have been widely used for epidemiological studies but were considered to have limited value as phylogenetic tools  <cit> .
very recently, however, the accuracy of spoligotyping for sublineage classification has been acknowledged  <cit>  and is <dig> has been shown to yield high phylogenetic resolution when used as target in a fluorescent amplified fragment length polymorphism technique  <cit> .
similarly, the is-seq strategy presented here resolves precise insertion site patterns characteristic of various strain families, thus providing robust phylogenetic information.
it remains to be seen, however, if the signatures identified herein will remain family-specific as more strains are analyzed from lineages not included in the present study.
although no correlation was found between is <dig> insertions present in different families and particular biological traits, we cannot rule out the possibility that clade specific insertions can alter important biological features of a strain and hence contribute to the spread of a new clade.
this study can serve as basis for additional analyses regarding the possible effects of these elements in clinical isolates and whether they may influence clonal spread by playing a role in strain fitness.
thus, this study presents a new method for efficiently typing naturally occurring insertion sequences in hundreds of isolates simultaneously.
it also expands the current inventory of in vivo interrupted loci, gives new insights into the genetic requirements in circulating strains and opens the way to understanding the consequences of insertions and the role of is <dig> in the biology and pathogenesis of m. tuberculosis.
furthermore, this approach can be used to study the natural history of insertion elements and transposons in the genomes of other pathogenic organisms where limited dna quantity and complex genome structure have made this pursuit challenging.
as costs of dna sequencing decrease, sequenced-based approaches will probably become more attractive for strain typing in a near future.
however, in a bacterial genome as stable as m. tuberculosis an approach that specifically targets variable regions, such as insertion sites that can provide sufficient phylogenetic information, may still prove more tractable than full genome re-sequencing of strains.
CONCLUSIONS
this article describes the simultaneous identification of is <dig> flanking regions in a large collection  of m. tuberculosis clinical isolates.
this is an innovative approach that combines an efficient barcoding strategy with a 96-well plate automated library preparation and massive parallel sequencing in a single illumina sequencing reaction.
our results showed about  <dig>  different insertions in the whole panel of analyzed isolates, which is more than has been identified up to date.
among the interrupted genes we identified some that were proposed to be essential for virulence and in vitro growth.
these results will be useful for researchers in the field of tuberculosis including the development of new drugs and vaccines.

 in spite of effective chemotherapy against tuberculosis, this disease is still a global health problem and a leading cause of death worldwide. the causative organism, mycobacterium tuberculosis, is an intracellular pathogen that has infected humans since ancient times. tuberculosis disease results from an intricate interaction between the host immune system’s efforts to control the infection and the pathogen’s ability to grow and persist within the host. thus, infection with the tubercle bacillus has variable outcomes that range from sterilizing immunity to active tb. active disease occurs in only 5– percent of immunocompetent individuals, with pulmonary tuberculosis being responsible for transmission in the community. in most cases, however, the infection is controlled by the host’s immune system and can lead to the establishment of a long-term latent infection which can reactivate later in life. there is no evidence of recent genetic exchange in m. tuberculosis which is thus considered to have a clonal population structure, with strains having almost identical nucleotide sequences. however, there is substantial intra-species genetic diversity that can affect disease outcome. the insertion element is is an important source of genomic variability for m. tuberculosis and is distributed in different positions across the bacterial chromosome. the copy number of this element is highly variable, but most strains contain between and instances. as a result, is has been widely used as an epidemiological marker in tuberculosis by dna fingerprinting using restriction fragment length polymorphism. these insertions also represent “in vivo” transposition events that provide information regarding genes required for human infection and disease. in addition, the persistence of this element in the m. tuberculosis genome raises the possibility that it might also drive phenotypic variability or affect strain fitness. it has been demonstrated that is transposition can be stimulated by specific stress conditions and thus contribute to genetic diversity in circulating m. tuberculosis clinical isolates. is insertions can affect gene expression by interrupting protein-coding genes, by mediating recombination events that result in deletions and inversions, and by up-regulating the expression of nearby genes due to a promoter located within the transposable element. thus, although is is primarily thought of as a valuable epidemiological tool, its prevalence and effect on genome function prompted us to take a deeper look at the distribution and patterns of is insertions by conducting a broad survey of circulating strains. a previous examination of is insertion sites in m. tuberculosis isolates, which was accomplished by pcr amplification of the target sequences followed by cloning and sequencing, identified unique insertions sites in only genes, suggesting that a large gene repertoire is required for human infection. however, the technical difficulties and biases due to the amplification and cloning associated with this approach prevented the analysis of a large number of isolates and the detection of all of the insertions present in each strain. to overcome these limitations and assess the extent of naturally occurring is insertions during infection and transmission, we used illumina® sequencing technology to identify is flanking regions in more than m. tuberculosis isolates from representative collections of clinical isolates from europe and south america. this high-throughput approach obviates the need for cloning and allows analysis of many strains in a cost-effective manner. we were able to identify approximately is insertion sites that interrupt almost genes, which together represent <10% of the m. tuberculosis genome. thus, we substantially increased the survey of naturally occurring in vivo is transpositions, providing new insights regarding their distribution, diversity and role in strain classification  insertion site sequencing we developed a method for the simultaneous identification of insertion sites in up to different strains by barcoding samples and using illumina® sequencing technology to sequence the genomic dna flanking both ends of an insertion sequence. we selected m. tuberculosis strains for analysis of the is insertion sites. this set included the reference strain m. tuberculosis h37rv and six other strains that were used to validate the method by verification of insertion sites using pcr. all of the strains included in this study were previously characterized by is rflp and spligotyping. *ppv: positive predicted value. for each reference strain the number of rflp observed bands is reported as well as the number of is sites identified using is-seq by one or both termini of the insertion and the number of confirmed sites by pcr. sensitivity and positive predicted value is calculated as indicated in methods. to generate a sequencing library we randomly sheared genomic dna and ligated barcoded adapters to the fragment ends. fragments that include the terminal end of the is were then amplified in a nested pcr using primers that incorporate additional barcodes. finally, the amplified products were pooled, purified and sequenced. the resulting fragments contain either the 3’ or 5’ terminal sequence of is the flanking genomic sequence, and two different barcodes at each terminus. since most of the enzymatic steps are “add only” reactions, there is no material loss due to purification steps and the protocol can be easily automated and performed in 96-well plates. furthermore, reactions are pooled after both barcodes have been added, so that the gel purification of hundreds of libraries can be performed using a single gel lane. using this protocol, we obtained more than million “paired-end” sequences using an illumina® ga ii machine. the sequence reads contained an even distribution of barcodes, showing that pooling before the gel extraction without individual quantification introduced no bias on the distribution. we were able to identify sequences for all of the strains with an average of reads per strain. to increase the precision of mapping, the short reads were sorted by sample and assembled to generate contigs of up to nucleotides before being mapped to the reference m. tuberculosis genome. all assembled contigs were mapped to at least one of the available genomes. the sequences that did not assemble in contigs were further mapped as single reads to the different genomes and less than 5% of these reads could not be mapped. based on an estimate of the average sequence coverage obtained per is element, a cutoff point was defined to remove strains with low coverage, which are likely to produce unreliable results. after this filtering, strains were retained for analysis. we included a negative control species in this study that lacks is and obtained only few and scattered reads without significant mapping. these results demonstrate the efficiency and specificity of the protocol and primer design. accuracy of the method the accuracy of is-seq was assessed using seven m. tuberculosis strains for which insertions had been previously characterized by rflp fingerprinting. is-seq identified these as well as additional insertions. of these additional insertions, all but two were confirmed by pcr. these two false positives corresponded to two insertions in the h37rv reference strain that were mapped to only one terminus of the is element. thus, by requiring sequence from both ends of the is element 100% specificity could be achieved, at the cost of sensitivity. the identification of additional insertion sites using is-seq was not altogether surprising given the limited resolution of the rflp fingerprinting method. based on the analysis of these strains, the is-seq method had a sensitivity of 100% and a positive predictive value of 98%. it was also interesting to find that while the col-h37rv strain used here had copies of is the sequenced genome of m. tuberculosis h37rv had is elements, with an extra copy being part of a doublet in positions, and,, highlighting differences in the genetic makeup of these reference strains. the reproducibility of the method was assessed by including two sets of duplicate samples of the same dna; the agreement on the number and insertion position within each pair of dnas was greater than 90%. a hallmark of the is element insertion is the generation of a 3- bp duplication at the insertion site. we expected to identify this duplication in those cases where we were able to map the insertion site of the is element on both termini. in % of these cases we identified a bp or bp duplication. taken together, these results demonstrate that the is-seq method accurately and reproducibly identifies and maps is insertion sites. genomic distribution of the is element the is-seq method identified insertion sites in the strains analyzed. these represented different insertions that affected loci, 66% of which were in coding regions and 34% were intergenic. this represents a significant bias towards non-coding regions, since intergenic regions account for only ~10% of the m. tuberculosis genome. an analysis of the distribution and frequency of independent insertion sites throughout the genome also showed a significant location bias. a poisson test for the distribution of insertion sites showed that there were forty-four kb-windows with significantly more insertions than expected. when kb windows were used, regions with significantly more or less insertions were also identified. an analysis of these regions revealed no significant correlation with gc content. we found multiple insertion sites in 40% of the affected loci. of these loci, contained or more independent sites and thus were labeled as “hotspots”, with the “hottest” locus containing different insertion elements. several of these “hotspots” have been previously described, such as mmpl, which is thought to be involved in fatty acid transport, or plcd, which encodes a phospholipase c, and intergenic regions like dnaa-dnan or the region between genes rv2813-. this latter locus corresponds to a cluster of regularly interspaced short palindromic repeats known as the direct repeat region, which is the template for the spoligotyping technique. we identified insertions in different sites within the dr region, one of which was present in of the evaluated strains, including the reference h37rv genome. interestingly, this is the only is present in the vast majority of mycobacterium bovis isolates, suggesting that it could have been present before many of the modern lineages diverged from each other. we therefore refer to this as an ancestral point of insertion. this insertion was absent in the beijing strains evaluated in our study, which is congruent with the loss of spacers 1– in the dr region that characterizes this m. tuberculosis family. to test if our analysis was close to saturation in terms of identifying all possible genes that could be naturally interrupted by the is element, a rarefaction analysis was carried out. the lack of a plateau indicating saturation in mutable sites shows that additional insertions in coding sequences could still be identified in circulating m. tuberculosis strains. this study was done with representatives of the major lineages circulating in europe and latin america and falls short of representatives of strains prevalent in africa and asia that may harbor distinct insertion sites. a regression analysis on the rarefaction curves estimated that between and genes could support insertion sites, consistent with previous estimates of the number of possible in vivo mutations. the presence of mutable sites results from an ongoing evolutionary process that in m. tuberculosis might not necessarily have reached saturation. thus, even if we sampled all available isolates there is no guarantee that we could observe such saturation of in vivo insertions. analysis of disrupted genes we next characterized the genes disrupted by is insertions by assigning them to functional categories based on the ncbi clusters of orthologous groups and the tuberculist database. disrupted genes were broadly distributed across functional categories, with a significant over-representation of insertions in genes belonging to three categories: 1) pe/ppe, 2) genes in cog from the category “unknown” and 3) cell wall/membrane/envelope biogenesis genes. the pe and ppe gene families encode glycine-rich proteins that contain repeats of the motifs pro-glu or pro-pro-glu and represent about 10% of the predicted proteins in m. tuberculosis. within the pe/ppe gene family, the majority of interrupted loci occurred in ppe genes. under-represented cog categories included genes involved in lipid and amino acid transport and metabolism, energy production and conversion, as well as the category defined as virulence, detoxification and adaptation in tuberculist. proteins belonging to the pe/ppe families were extracted to a separate category, as they constitute an important family of proteins in m. tuberculosis. the genes for the transposase is were removed from the category of replication, recombination and repair. cumulative length of all the genes in a given category. the probability of under or over-representation of a given functional category is dependent on both the number of genes and their length. represents the number of independent insertion events identified in genes of a given category. probability of over or under-representation of insertion sequences interrupting genes of a given category. in bold, categories with significant over-, under-representation after bonferroni correction. bonferroni corrected threshold = e- this category is not part of cog; it is defined in tuberculist. we next compared the set of orfs harboring an is element with previous studies that identified m. tuberculosis genes required for growth in vitro or in vivo. twenty-seven genes previously identified as being essential under different experimental conditions were found to contain is elements. seven of these genes contained the insertion element in the terminal 20% of the gene, which although less likely to affect function could still do so in proteins with c-terminal active sites or domains. to rule out the possibility that these observations were due to an artifact of the is-seq method, we validated these insertion sites by pcr. twenty-four of the insertions were confirmed, and one locus failed to amplify, probably due to the limited amount of dna used as starting material. rv3473c indicates the position of the insertion within the orf as the percentage of gene upstream of the insertion. 2) insertions verified using pcr: +, pcr positive; -, pcr negative. is family signatures to see how well these is patterns correlated with known strain classifications, we examined the relationship between is distribution and m. tuberculosis strain families. the use of specific is insertions as markers for different m. tuberculosis lineages has been described. for example, the beijing lineages have been described to contain an is insertion in the intergenic region between dnaadnan and the spoligotype characteristic of the haarlem lineage is partially due to a characteristic is in the dr region. despite the limitation of using the is to accurately discriminate between the different m. tuberculosis phylogenetic lineages, we propose that the precise locations of is insertions could be useful as markers within a lineage and, furthermore, could help to identify previously uncharacterized clades within families. we therefore used a mutual information metric to identify the co-occurrence or absence of insertion sites in order to generate a dendrogram using of our strains, which belonged to the m. tuberculosis euro-american lineage and the beijing family, which was also used as an outgroup. we also included the publicly available genomes that belong to either lineage or. the tree generated distinguished the main lineage families, and subdivisions within families, such as the h h and h haarlem sub-families defined by spoligotyping. for the lam family, two main clusters were formed, one composed of lam lam lam and the second of lam lam lam and tuscany. the t family also showed sub-clusters that did not reflect the t sub-families described by spoligotyping, but were consistent with findings using snps for classification that suggested that t sub-families are still not well defined. by analyzing the patterns of insertion events in the different families, we identified family-specific insertions, consisting of few insertions per family. these signature is sites helped to understand inconsistencies in spoligotyping. in particular, strains ambiguously classified by spoligotyping as either ‘haarlem orphan, h and h and only as haarlem by rflp, were found to contain a h subfamily characteristic insertion site and a h subfamily spoligotype pattern in the dr region, which may have occurred independently as a product of recombination between two is interestingly, we were able to identify the insertions in the dr region that are responsible for the spoligotype patterns of the different members of the haarlem family. two haarlem strains were found to be ural strains, consistent with recent reports, and the tuscany subfamily was closely related to the lam subfamily rather than the t family, as had been recently suggested. also, the previously described insertion in the dnaa -dnan intergenic region in beijing isolates was identified in all strains. when these family-specific insertions were used to classify the strains used in this study, family or subfamily information was provided for strains compared to that were classified by spoligotyping. these signatures were also used to classify of the available complete m. tuberculosis genomes. some strains were assigned to specific families, as was the case of the reference h37rv and the h37ra strains that were assigned to the t family, and the f and haarlem strains that were identified as belonging to lam and haarlem respectively. other strains were classified as lam, lam, x or beijing; all in 100% agreement with previously suggested families, when available. the inability to classify the remaining cases is due to the fact that they belong to lineages not significantly represented in our study. these results show that classification of tuberculosis strains by the use of signature is sites is consistent with snp classification and therefore highly reliable. to address the possibility that insertions could appear independently in different lineages because they confer some selective advantage, such as drug resistance, we performed a conditional mutual information analysis to see if any insertion sites correlated significantly with the epidemiological data collected for each strain. no association was found between insertion sites present in different phylogenetic families and a specific biological trait, suggesting that it is unlikely that exact, yet independent insertions could occur in different families that confer a fitness advantage. furthermore, in all cases where a specific insertion sequence was present in multiple strains, the strains containing that insertion belonged to the same phylogenetic lineage, suggesting that occurrence of independent insertions in the same nucleotide position but in different lineages by convergent evolution, although plausible, is very unlikely. instead, hotspots of insertions occurring in a single locus within few nucleotides of distance were observed. discussion the is-seq method was effective at detecting, in a single reaction, multiple is insertions sites in over m. tuberculosis isolates from diverse geographical origins and phylogenetic lineages. this was achieved using a novel approach with the following features: i) the use of two molecular barcodes per sample to multiplex samples on an illumina® sequencer, allowing the analysis of a combinatorial number of samples with just different barcodes; ii) a high-throughput library construction protocol that uses minimal amounts of starting material, and iii) a novel computational strategy that maximizes the efficiency of sequence identification by assembling reads prior to mapping against the m. tuberculosis sequenced genomes. this study, combined with the advances on sequencing technologies, provides a highly cost-efficient method for m. tuberculosis strain identification and typing that could be extended to other target amplifications and sequencing. while the use of bp barcodes allows multiplexing of up to samples per sequencing reaction, bp barcodes may allow pooling up to samples. additionally, the requirement of just reads per insertion site for accurately mapping, combined with new technologies such as miseq that achieve up to million reads in just a matter of hours, can revolutionize the field of tuberculosis strain characterization. is-seq effectively detected loci with multiple independent insertions or insertions in the same location in many strains, and identified new and previously characterized hotspots and insertions that could be considered ancestral. although an is element present in the same locus could be due to convergent evolution, which in time could be reflected as insertion hot spots, our observations suggest that insertions at the same nucleotide position are more likely due to a common ancestral insertion. strategies similar to is-seq have been recently developed for identifying microbial transposon insertions sites generated during in vitro mutagenesis and endogenous transposon elements in the human genome. however, our method differs because it can identify multiple insertions sites per strain and simultaneously analyze hundreds of different strains per reaction. the advantage of is-seq results from a more streamlined library preparation protocol and a more efficient barcoding strategy. the insertion sites identified in these clinical m. tuberculosis isolates significantly extends previous surveys by identifying more that new sites and provides, for the first time, a comprehensive overview of the gene repertoire that can be interrupted in vivo by is rarefaction curves indicated that additional genes could still be identified, especially in strains that were not significantly represented in our dataset. however, the number of loci predicted to sustain in vivo insertions was close to the maximum of about genes previously predicted based on an analysis of strains. it is also consistent with the observation that the number of naturally occurring in vivo mutations is much lower than the genes predicted to be non-essential using transposon mutagenesis approaches, probably due to the detrimental effect of cumulative insertions on in vivo fitness. the smaller gene set that can be inactivated in vivo when compared to in vitro studies, tied to identification of insertions in previously predicted essential genes, indicates evident differences in the conditions imposed on circulating strains as they adapt to the human host. in this respect, the overrepresentation of insertions in pe/ppe, cell wall and membrane biogenesis genes and orfs of “unknown function” could reflect adaptation as bacteria alter envelope components in response to pressure imposed by the host environment. in particular, the presence of is elements in the pe/ppe gene family reinforces their previously proposed role in immune evasion and is consistent with previous reports. however, the over-representation of certain functional categories could also be due to redundancy in function that might alleviate the effect of gene disruption. in contrast, underrepresented gene categories, such as energy production and conversion, highlight their importance in processes required for bacterial virulence and survival. finally, the identification of is insertions in genes reported to be essential is not implausible, since insertions as well as deletions of some of these genes have already been described. in addition to gene disruption, the is insertions identified here can also result in transcription up-regulation of downstream genes, as has been reported for an insertion in the promoter region of the two-component system gene phop. a similar insertion in the phop promoter was observed here in a fast evolving clade of haarlem in latin america. this survey therefore broadens the list and emphasizes again the difference in requirements between strains grown and studied under laboratory conditions and those subjected to the pressure of in vivo passage and processes of natural infection in the human host. interestingly, the strain classification based on the distribution of is using is-seq was congruent with that based on spoligotyping and snps for sequenced genomes. both is rflp and spligotyping have been widely used for epidemiological studies but were considered to have limited value as phylogenetic tools. very recently, however, the accuracy of spoligotyping for sublineage classification has been acknowledged and is has been shown to yield high phylogenetic resolution when used as target in a fluorescent amplified fragment length polymorphism technique. similarly, the is-seq strategy presented here resolves precise insertion site patterns characteristic of various strain families, thus providing robust phylogenetic information. it remains to be seen, however, if the signatures identified herein will remain family-specific as more strains are analyzed from lineages not included in the present study. although no correlation was found between is insertions present in different families and particular biological traits, we cannot rule out the possibility that clade specific insertions can alter important biological features of a strain and hence contribute to the spread of a new clade. this study can serve as basis for additional analyses regarding the possible effects of these elements in clinical isolates and whether they may influence clonal spread by playing a role in strain fitness. thus, this study presents a new method for efficiently typing naturally occurring insertion sequences in hundreds of isolates simultaneously. it also expands the current inventory of in vivo interrupted loci, gives new insights into the genetic requirements in circulating strains and opens the way to understanding the consequences of insertions and the role of is in the biology and pathogenesis of m. tuberculosis. furthermore, this approach can be used to study the natural history of insertion elements and transposons in the genomes of other pathogenic organisms where limited dna quantity and complex genome structure have made this pursuit challenging. as costs of dna sequencing decrease, sequenced-based approaches will probably become more attractive for strain typing in a near future. however, in a bacterial genome as stable as m. tuberculosis an approach that specifically targets variable regions, such as insertion sites that can provide sufficient phylogenetic information, may still prove more tractable than full genome re-sequencing of strains.  this article describes the simultaneous identification of is flanking regions in a large collection of m. tuberculosis clinical isolates. this is an innovative approach that combines an efficient barcoding strategy with a 96-well plate automated library preparation and massive parallel sequencing in a single illumina sequencing reaction. our results showed about different insertions in the whole panel of analyzed isolates, which is more than has been identified up to date. among the interrupted genes we identified some that were proposed to be essential for virulence and in vitro growth. these results will be useful for researchers in the field of tuberculosis including the development of new drugs and vaccines. 
although is has been used extensively as an epidemiological marker, the identification of the precise chromosomal insertion sites has been limited by technical challenges. strains were classified into families using insertion sites, and high agreement with previous studies was observed. the insertions mapped to both coding and non-coding regions, and in some cases interrupted genes thought to be essential for virulence or in vitro growth. the insertion element is is one of the main sources of genomic variability in mycobacterium tuberculosis, the etiological agent of human tuberculosis. this high-throughput is-seq method, which can also be used to map insertions in other organisms, extends previous surveys of in vivo interrupted loci and provides a baseline for probing the consequences of disruptions in m. tuberculosis strains. when validated using reference strains, the method had 100% specificity and 98% positive predictive value.
 the insertion element is is one of the main sources of genomic variability in mycobacterium tuberculosis, the etiological agent of human tuberculosis. although is has been used extensively as an epidemiological marker, the identification of the precise chromosomal insertion sites has been limited by technical challenges. here, we present is-seq, a novel method that combines high-throughput sequencing using illumina technology with efficient combinatorial sample multiplexing to simultaneously probe clinical isolates, identifying almost all the flanking regions of the element in a single experiment.  we identified a total of is flanking regions on the different isolates. when validated using reference strains, the method had 100% specificity and 98% positive predictive value. the insertions mapped to both coding and non-coding regions, and in some cases interrupted genes thought to be essential for virulence or in vitro growth. strains were classified into families using insertion sites, and high agreement with previous studies was observed.  this high-throughput is-seq method, which can also be used to map insertions in other organisms, extends previous surveys of in vivo interrupted loci and provides a baseline for probing the consequences of disruptions in m. tuberculosis strains. 