Spoligotyping analysis of Mycobacterium tuberculosis in Khyber Pakhtunkhwa area, Pakistan

Introduction

Tuberculosis (TB) has become a major public health problem in Pakistan over the past few decades, and was declared a global emergency by the WHO in 1993. Various national as well as international authorities are highly concerned and are contributing their efforts to abolish TB. Pakistan is among 22 countries identified by the WHO for its high burden of TB cases reported annually. Drug resistance is a major obstacle to TB control in Pakistan.^1,2 Millions of Afghan refugees have seasonal traffic between the two countries, which suggests an epidemiological linkage and evolutionary genomic characteristics of Mycobacterium tuberculosis (MTB) strains in Pakistan. The genome of MTB contains variable copies of 36 bp direct repeats (DRs) in a specific locus, the CRISPR (clustered regularly interspaced short palindromic repeats) locus. These DR loci are interspersed with 34–41 bp conserved sequences called spacers. More than 100 spacers have been identified so far; however, only 43 spacers are used for genotyping while the others contribute only a small amount of information.³ This set of 43 unique spacer sequences has been used to classify MTB strains based on the pattern of the presence or absence of DRs.⁴ Variation in the number of spacers contributes to combinational or individual effects of the genomic process, such as homologous recombination, DNA replication or transposition of the IS6110 element in the DR locus. The function of CRISPR in bacterial metabolism and its relationship with the spacers have not yet been well defined in MTB;⁵ however, it is evident from the literature that certain patterns of DRs are linked to evolutionary processes.³

Spoligotyping is one of the most highly adopted PCR-based methods for genotyping analysis and is considered both reliable and cost-effective as a first-line genotyping method. Various new techniques have been proposed for spoligotyping analysis, including the DNA microarray,⁶ PixSysn QUAD 4500 Microarrayer,⁷ TB-SPRINT⁸ and hydrogel microarray.⁹ TB-SPRINT is a microbead assay that provides a spoligotyping profile and simultaneous detection of common mutations in the hot-spot regions of rpoB, katG and inhA genes to evaluate drug resistance to rifampicin and isoniazid in MTB. TB-SPRINT is a useful tool for epidemiological studies of TB in resource-limited countries because of its high throughput and acceptable cost.⁸ The aim of the current study was to investigate the epidemiological link between MTB strains and their genotypic families found in Pakistan using the spoligotyping method.

Methods

This study was approved by the research and ethics committee of the Provincial TB Control Program Khyber Pakhtunkhwa (TB/KP/R&D-25-17). Samples for the current study were collected at the Provincial TB Reference Laboratory, Peshawar, which is equipped with a biosafety level-III (BSL-III) facility. Samples were digested and decontaminated using a standard N-acetyl-L-cysteine sodium hydroxide (NALC-NaOH) method. Decontaminated samples were inoculated on Lowenstein–Jensen (LJ) medium and a mycobacterium growth indicator tube (MGIT). Positive growth in the tubes was confirmed using a TBc Identification Test device (ref. 245159; Becton Dickinson, Franklin Lakes, NJ, USA). All confirmed mycobacterial isolates were subsequently processed for genotyping analysis targeting the DR locus in MTB, as described by Kamerbeek et al.⁴ DNA was extracted from the culture using the cetyl–trimethyl-ammonium-bromide (CTAB) method.⁴

Spoligotyping assay was performed on all samples using a Luminex 200 (Luminex Corp., Austin, TX, USA) instrument.⁸ Spoligotype patterns were entered in binary format into the SITVIT2, a proprietary database of the Institut Pasteur de la Guadeloupe, which is an updated version of the SpolDB4 and SITVIT-WEB databases.^10,11 Clinical isolates were initially grouped in one of nine lineage and sublineages: L1/EAI, L2/Beijing, L3/CAS, L4.1.1/X, L4.1.2/Haarlem, L4.3/LAM, L4.4/S, L4/T-others, L5 and L6 (Mycobacterium africanum WA1 and WAS2).¹¹ The selected clinical isolates that matched with a published spoligo-international type (SIT) were assigned with the corresponding SIT label, while strains with an unknown pattern were declared “orphans”. Lineages were deduced based on the latest next-generation sequencing-based taxonomy when correspondence between whole-genome-based and spoligotyping-based taxonomy was possible.^12,13 All strains whose pattern did not fall within the above-mentioned lineages were declared new patterns. Furthermore, the patterns of all strains originating from Pakistan were collected from the SITVIT2 database. Evolutionary relationships between the unclustered spoligotype patterns and strains reported from Pakistan were determined using MEGA7 software.¹⁴

The genetic diversity was laid out as a maximum parsimony (MP) tree construct using default settings in MEGA7. The phylogenetic tree was constructed using the subtree–pruning–regrafting (SPR) algorithm, taking the 43 DR pattern as a sequence.¹⁴

Results

In total, 200 isolates of MTB were analyzed by the spoligotyping method. Of these isolates, 166 produced valid and interpretable results. Patterns of 145 isolates were found in the international database (SITVIT2) and 21 unique patterns (13%) were declared new types (Table 1). The results showed that lineage 3/Central Asian strain (L3/CAS) was the predominant genotype lineage with 72.8% (n=121), followed by L2/Beijing with 5.4% (n=9) and L4/T with 4.2% (n=7), and 4.8% (n=8) were other unknown types (Figure 1).

Table 1 Spoligotypes detected in Mycobacterium tuberculosis isolates from KP, Pakistan

Figure 1 Graphical analysis of 166 Mycobacterium tuberculosis isolates collected from different areas of Pakistan. Note: Lineage 3/Central Asian strain (L3/CAS), shown in black, is the dominant spoligotype family. New unclustered strains are shown in the black dotted segment. Less frequent spoligotypes are grouped in the “Others” category in the gray segment.

Cluster analysis of widely distributed lineages

Detailed analysis of L3/CAS resulted in 11 different clusters or sub-families. Among the sub-families, L3/CAS1-Delhi (SIT26) was the predominant genotype (n=78), followed by L3/SIT1401 (n=5) and L3/SIT794 (n=4). The pattern of one member of the L3/CAS family that did not find a match in the SITVIT2 database was declared as L3/CAS-Orphan. All strains belonging to the L2/Beijing family were grouped in two clusters, SIT1 (n=8) and SIT941 (n=1). Members of L4/T-other spoligotypes were grouped in five different clusters (SIT51, SIT53, SIT273, SIT732 and SIT189) and one other T orphan. It was also found that EAI1-SOM, EAI5, EAI6-BGD1, H1, LAM9-ZWE, Manu, EAI3-IND and X1 were rarely distributed in this region.

All 21 strains with unique spoligotyping patterns were grouped into 13 clusters based on the similarity index of the MP test. The new clusters were further analyzed using different methods to determine any potential evolutionary relationship to previously published genotypes. The results indicated that the majority of the unique patterns were sublineages of the L3/CAS family, except for some orphans (KP1–KP3) that were more closely related to L4. The detailed phylogenetic tree is represented in Figure 2. These spoligotyping patterns were unique not only to Pakistan but also to the international database. Therefore, no SIT has been assigned so far (Table 1).

Figure 2 Maximum parsimony analysis of new unclustered spoligotypes.

The tree shows the relationship of new unclustered spoligotypes in this study to the previously published spoligotype patterns from Pakistan.^15–19 The new spoligotype patterns in this study were numbered Orphan KP1 to Orphan KP13.

Discussion

Many techniques have been used for MTB genotyping, including restriction fragment length polymorphism (RFLP), amplified fragment length polymorphism (AFLP), mycobacterial interspersed repetitive unit–variable number tandem repeat (MIRU-VNTR) and spoligotyping. The spoligotyping technique has been used worldwide for genotyping and as a first-line method in epidemiological studies of MTB, in conjunction with VNTR analysis. It was originally developed by Kamerbeek et al, based on identical DRs interspersed by short spacers of variable length.⁴ A public database, SITVIT2, is available for the analysis of spoligotyping patterns, which consists of 111,635 clinical isolates from 169 countries of patient origin.²⁰ The distribution of these spoligotypes was analyzed among the 20 studied subregions. It is evident from the literature that these genotypes are specific to certain geographical regions; for example, L4.3/LAM is geographically widely distributed in South America, followed by L4.1/T and L4.1.2/Haarlem lineage, at 49.3%, 26.7% and 15.7%, respectively. Similarly, L1/EAI is geographically distributed in South-East Asia and L3/CAS in Middle-East Asia.^11,21

Our study represents the first molecular characterization of TB bacilli in the study area and may be useful for better management of TB and drug resistance.²² Considering the spoligotyping patterns of the majority of the strains, the current study confirmed a diverse population structure with a total of 55 different types, whether clustered or not. However, by far the most dominant lineage was L3/CAS-1 Delhi, which belongs to the L3/CAS lineage. Overall, the distribution of genotypes in Khyber Pakhtunkhwa was in line with other spoligotyping analysis in Pakistan.^16,19,21 However, we do not report any genotype of L4.5 “NEW-1”, an emerging genotype in Pakistan, since no SIT127 types were observed.²³ Spoligotyping, in some cases, is not able to identify L4.5, which is better characterized by RD122, specific VNTR signatures and specific single-nucleotide polymorphisms.^13,23 A rather important population of 12% of strains could not be matched to the internationally defined genotypes in the database and was declared “new”. In previous investigations, the spoligotyping patterns of new genotypes were underestimated and left unstudied.^16,19 In-depth analysis revealed clear identification of L3/CAS by the presence of the first three spacers followed by a set of four absent spacers (4–7) and two more spacers (33–34). We confirmed that the majority of the new spoligotypes in this study were likely to have evolved from L3/CAS, characterized by the absence of spacers 4–7 and 23–34, as changes in the DR locus are most likely to occur via consecutive irreversible deletions of either single or contiguous blocks.^3,24

The current study provides basic information on the geographical distribution of MTB strains in Khyber Pakhtunkhwa, Pakistan. It has some limitations for the time being, since spoligotyping-based clusters were not further studied using VNTR. However, such comparisons of spoligotypes may not be representative of the true phylogenetic relationships owing to genetic convergence. By contrast, DNA sequencing can be used to determine the true evolution.

Conclusion

We found that L3/CAS is the predominant lineage, which is likely to be epidemiologically linked in many cases, and the majority of the unclustered L3 strains are likely to have evolved from pre-existing L3/CAS strains based on the MP analysis tree construct. Other tools, such as MIRU-VNTR in combination with spoligotyping, would help us to further explore the molecular epidemiology of MTB in our study area.