Whole Genome Analysis of Extensively Drug Resistant Mycobacterium Tuberculosis Strains in Peru&nbsp;

doi:10.21203/rs.3.rs-122146/v1

Download PDF

Research Article

Whole Genome Analysis of Extensively Drug Resistant Mycobacterium Tuberculosis Strains in Peru

https://doi.org/10.21203/rs.3.rs-122146/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Peru has the highest burden of multidrug-resistant tuberculosis in the Americas region, and since 1999, the annual number of extensively drug-resistant TB (XDR-TB) cases has been increasing, becoming a public health challenge. The objective of this study was to perform genomic characterization of Mycobacterium tuberculosis strains obtained from peruvian patients with XDR-TB diagnosed from 2011 to 2015 in Peru.

Whole genome sequencing (WGS) was performed on 68 XDR-TB strains from different regions of Peru. 58 (85.3%) strains came from the most populated districts of Lima and Callao. Concerning the lineages, 62 (91.2%) strains belonged to the Euro-American Lineage, while the remaining 6 (8.8%) strains belonged to the East-Asian Lineage. Most strains (90%) had high-confidence resistance mutations according to pre-established WHO-confident grading system. Discordant results between microbiological and molecular methodologies were due to mutations outside the hot spot regions analyzed by commercial molecular assays (rpoB I491F and inhA S94A). Cluster analysis using a cutoff ≤ 10 SNPs revealed that only 23 (33.8%) strains were recent transmission events.

This study highlights the relevance and utility of WGS as a high-resolution approach to predict drug resistance, analyse transmission of strains between groups, and determine evolutionary patterns of circulating XDR-TB strains in the country.

Epigenetics & Genomics

Peru

extensively drug-resistant TB (XDR-TB)

Lima and Callao

WHO-confident grading system

Tuberculosis (TB) is a preventable and curable disease that is one of the top 10 causes of death in the world ¹. Resistance to drugs used in the treatment of TB is a major threat to the strategies that are being deployed to control and eliminate the multidrug-resistant TB (MDR-TB) and extensively drug-resistant TB (XDR-TB) ^1–3. Both are severe forms of this pathology that are becoming more prevalent through years and require an expensive and prolonged treatment that produces greater morbidity, toxicity and mortality ^1,4. In 2018, worldwide it was reported that 3.4% of new TB cases and 18% of previously treated cases were diagnosed with MDR-TB, of which 6.2% additionally showed to present XDR-TB, which has the lowest treatment success rate (39%) compared to other forms of TB ¹. XDR-TB is present in the region of the Americas among which Peru is the country with the highest burden ^1,5. The first Peruvian case of XDR-TB was detected in 1999 ⁶, and throughout the years the cases had progressively increased reporting a cumulative number of 944 cases until 2017 ^7,8.

Drug-resistance of mycobacterial strains is detected through genotypic or phenotypic laboratory tests that detect the presence of DNA mutations conferring resistance or the growth of Mycobacterium tuberculosis (MTB) in the presence of antituberculous drugs, respectively. However, these methodologies are limited to the analysis of a limited number of resistance genes or to the slow mybacterial duplication time, making it difficult to obtain complete and rapid results. In addition, for relatedness determination, analysis of the polymorphic length restriction fragments of the IS6110 gene (IS6110 RFLP), analysis of spacer oligonucleotide typing (Spoligotyping) and analysis of variable number of tandem repeats of interspaced repetitive mycobacterial units (MIRU-VNTR) were implemented globally. These methodologies have also been applied in Peru for exploration of the genetic diversity of drug-resistant TB strains ^9–11. However, the problem using these conventional techniques for genotyping is that they explore polymorphic genetic regions that cover less than 1% of the mycobacterial genome, limiting their power of differentiation between different strains ^12,13.

Recently, the revolution of Next Generation Sequencing (NGS) technology and its wider availability have allowed the use of Whole Genome Sequencing (WGS) to provide information about speciation, drug resistance prediction and better determination of relatedness for epidemiologic purposes ^14–16. In this way it is possible to obtain a greater amount of information that allows a complete characterization and discrimination of strains with repeated or ambiguous conventional genotypic patterns ^17,18.

High genomic resolution through WGS has not been performed in the Peruvian TB-XDR strains. Therefore, the aim of the present study is to characterize the genetic variability of the XDR-TB strains circulating in Peru, through the NGS-based WGS analysis. Likewise, we performed a general approximation of their molecular epidemiology and determined the evolutionary relationships of the XDR-TB Peruvians strains.

Patient characteristics

68 strains were included. Five (7.4%) strains were obtained from patients diagnosed in 2011, two (2.9%) in 2012, 18 (26.5%) in 2013, 23 (33.8%) in 2014 and 20 (29.4%) in 2015. Likewise, 38 (55.9%) strains were initially obtained from men and 30 (44.1%) from women. Their age ranged from 15 to 78 years with a mean of 36 years (interquartile range 25.3–44.5). 51 (75.0%) cases had received a previous antiTB treatment (Table 1). 58 (85.3%) strains belonged to patients from Lima and Callao (51 from Lima and 07 from Callao). Regarding the strains of Lima, 30 (59%) came from the east zone (Lima este), 16 (31%) from the center and north (Lima ciudad), 3 (6%) from the south (Lima sur), and 2 (4%) from the other provinces (Lima provincias). In addition, Piura, La Libertad, Loreto, Ancash, Ucayali, Arequipa and Madre de Dios regions were represented by one strain each, whereas Ica region with two strains (Supplementary Table S2).

Table 1

**Characteristics of patients with XDR-TB**.
Characteristics	n [%]
Sex
Male	38 [56]
Female	30 [44]
Treatment History
Previously treated	51 [75]
Newly diagnosed	15 [22]
Unknown	2 [3]
Age
≤ 15 year.	1 [1]
16–30 year.	24 [35]
31–45 year.	25 [37]
> 45 year.	16 [24]
Unknown	2 [3]
Mean (IQR)	36.3 (25–45)
Second Line Inyectables*
Kanamycin resistant	11 [16]
Capreomycin resistant	14 [21]
Kanamycin & Capreomycin resistant	43 [63]

*All strains (n = 68) were additionally resistant to rifampicin, isoniazid and levofloxacin drugs.

IQR

interquartile range.

Sequencing and genome assembly

An average of 935,183 raw sequencing reads per fastq file was obtained. Two fastq files (forward and reverse) were generated for every sample. The minimum, maximum and average depth of sequencing obtained were 53x, 153x and 88x, respectively (Table S2). All strains had reads covering more than 99% of the H37Rv genome.

Antimicrobial resistance

All strains showed simultaneous phenotypic resistance to rifampicin, isoniazid and levofloxacin. However, they showed differences in second line injectable drugs resistance, as described in Table 1. Discordant results between phenotypic and genotypic methods were found for isoniazid (strain ‘XDR_28’) and rifampicin (strain ‘XDR_19’) showing only resistance outputs through phenotypic method (Table S1) and analysed elsewhere ¹⁹. Concerning rifampicin resistance, all strains had resistant mutations located at rpoB gene. From this, 67 (98.5%) were considered high-confident mutations located inside the rifampicin-resistance-determining region (RRDR) of the rpoB gene, while only one strain (1.5%) presented a mutation outside this region. The most frequent mutations were S450L, D435V and H445R. Only one strain had the double mutation H445N + S431R. The isoniazid-resistant phenotypes showed mutations in katG and inhA genes. There was one mutation in katG (S315T), whereas three differente mutations located in the promoter region (g-17t, c-15t and t-8c) and one mutation located in the codificant region (S94A) of inhA were found. There were four strains with double mutation, S315T + c-15t, and one with S315T + t-8c.

Levofloxacine resistance was predominately caused by mutations occurring in the quinolone-resistance-determining region (QRDR) of gyrA gene and contained nine different mutations, whereas gyrB showed showed only three different mutations. In the gyrA gene, the codon 94 showed the greatest variability (five different mutations). Only two strains presented mutations in both genes. Finally, one strain (1.5%) strain showed no mutations in any of these genes.

Resistance to kanamycin and capreomycin was driven by mutations occurring at rrs (a1401g, c1402t and g1484t), tlyA (all frame shifts) and and eis (c-14t) genes. However, no mutations were detected in two strains for the screened genes. Strains with exclusive resistance to kanamicine only showed the a1401g mutation of rrs gene, whereas two strains had no mutation in any of the three genes analyzed. Exclusive resistance to capreomycin was caused by several frameshifts’ mutations occurring in tlyA gene. However, there were three strains with no detected mutations in any of the analyzed genes. In general, 96, 85 and 90% (average 90%) of strains had high-confident mutations for resistance to rifampicin, isoniazid and both second line drugs, respectively (Table 2).

Several synonymous and nonsynonymous mutations located in additional resistant-associated genes were evidenced to be present together with the mutations described above (Supplementary Table S3).

Table 2

**Resistance-associated mutations in peruvian XDR MTB strains**.
Drug	Gene	Mutation	Confidence grading **	n (%)
Rifampicin	rpoB	Q432P	High	1 (1.5)
		D435V	High	22 (32.4)
		H445N, S431R	Minimal	1 (1.5)
		H445R	High	4 (5.9)
		H445S	No data	1 (1.5)
		S450L	High	38 (55.9)
		I491F	Minimal	1 (1.5)
	Total			68 (100)
Isoniazid	katG	S315T*	High	53 (77.9)
	inhA	c-15t	Moderate	7 (10.3)
		g-17t	No data	2 (2.9)
		S94A	No data	1 (1.5)
	katG + inhA	S315T* + c-15t	High	4 (5.9)
	katG + inhA	S315T* + t-8c	High	1 (1.5)
	Total			68 (100)
Levofloxacin	gyrA	G88C	High	2 (2.9)
		A90V	High	13 (19.1)
		S91P	High	2 (2.9)
		D94A	High	2 (2.9)
		D94G	High	32 (47.1)
		D94H	No data	3 (4.4)
		D94N	High	5 (7.4)
		D94Y	High	1 (1.5)
		A90V, D94H	High	1 (1.5)
		D94N, D94G	High	1 (1.5)
		D94N, D94A	High	1 (1.5)
	gyrB	S447F	No data	2 (2.9)
	gyrA + gyrB	G88A + T500P	No data	1 (1.5)
	gyrA + gyrB	D94A + E501D	High	1 (1.5)
	Mutation not detected			1 (1.5)
	Total			68 (100)
Kanamycin & Capreomycin	rrs	a1401g	High	34 (50.0)
		c1402t	High	1 (1.5)
		g1484t	High	1 (1.5)
	tlyA	S92fs	Generally high	2 (2.9)
		A111fs	Generally high	1 (1.5)
		H68fs, S92fs, Q202fs	Generally high	1 (1.5)
	rrs, eis, tlyA	A1401G, c-14t, R60fs	High	1 (1.5)
	Mutation not detected			2 (2.9)
Kanamycin	rrs	a1401g	High	9 (13.2)
Kanamycin	Mutation not detected			2 (2.9)
Capreomycin	tlyA	V198fs	Generally high	4 (5.9)
		G232D	Generally high	2 (2.9)
		L209dup	Generally high	1 (1.5)
		C86fs	Generally high	1 (1.5)
		S252fs	Generally high	1 (1.5)
		R133fs, L209dup	Generally high	1 (1.5)
		L139fs, V198fs	Generally high	1 (1.5)
	Mutation not detected			3 (4.4)
	Total			68 (100)

Codon numbering systems are according H37Rv genome (NC_000962.3). fs: frame shift, dup: duplication.

* AGC ➡ ACC

** According to WHO-NGS Technical guide ²⁰

Lineages and evolutionary relationships

Lineages analysis was performed in the 68 strains determining that 62 (91.2%) XDR strains belong to the Euro-American Lineage (Lineage 4) and 6 (8.8%) to the East-Asian Lineage (Lineage 2). 59 strains from Lineage 4 were able of sublineages classification, whereas all Lineage 2 strains were only represented by one sublineage. Moreover, in silico Spoligotypes analysis showed that within Lineage 4 there was a statistically higher frequency (p = 0.23) of the H, LAM and T families followed by family X. Likewise, the existence of only five (7.4%) strains with orphan spoligotypes and of three (4.4%) strains of unknown spoligotypes was evidenced (Fig. 1, Supplementary Fig. S1).

In general, the clustering of strains according to the patterns of spoligotypes agreed with the lineages and sublineages obtained; however, the strains labelled as belonging to the family T1 were assigned to different sublineages of Lineage 4. The northern strains of the country belonged to sublineages 4.1.1 (Piura, n = 1), 4.1.1.3 (La Libertad, n = 1) and 4.3.3 (Loreto, n = 1). In the center, the city of Lima contained 7 sublineages (except 4.1.1.3); likewise, the strains from the constitutional province of Callao were aligned in the sub-clustering groupings 4.3.3 (n = 3), 4.3.4.2 (n = 2), 4.1.2.1 (n = 1) and 2.2.1 (n = 1), evidencing a great diversity of strains that circulate in these areas. In the rest of the cities, strains of the subfamilies 4.3.3 (Ica and Ucayali, n = 2), 4.1.1 (Ancash, n = 1) and 4.3.4.2 (Huancavelica, n = 1) were found. Finally, in the South of the country, strains were found to belong to the sublineages 4.1.1.3 (Arequipa, n = 1) and 4.1.2.1 (Madre de Dios, n = 1) (Fig. 1).

The maximum likelihood phylogenomic tree confirmed the lineage classification and showed adittional subclassification for strains belonging to the same sublineage. Interestingly, there was no association between phylogenetic clades and specific resistant mutations conferring resistance. Finally, spoligotypes with SITs ‘Orphans’ and ‘Unknown’ could be characterized by evolutionary and lineage analysis. All strains with ‘Unknown’ SITs aligned within the sublineage 4.3.3 group, showing an evolutionary similarity with members of the H3 family, while the strains with ‘Orphans’ SITs were located in the groups belonging to sublineages 4.3.3 (n = 4), 4.1 (n = 1) and 4.1.2.1 (n = 1). A close evolutionary relationship was found between strains from Arequipa (‘XDR_10’) and La Libertad (‘XDR_05’), being the only representatives of the sublineage 4.1.1.3 in the entire country and exhibiting a greater degree of genomic differentiation with the other strains (Fig. 2). Global evolutionary relatedness showed that most of Peruvian XDR strains were grouped in defined clusters and had a close relatedness with XDR strains of Lineages 4 and 2 of european countries (Supplementary Fig. S2).

Transmission cluster determination

The analysis of pairwise genetic differences showed a high number of strains that differed by large amounts of SNPs. This genetic distance between strains varied from 5 to 1272 nucleotide differences with an interquartile range of 373 and a median of 772 (Supplementary Fig. S3). The WGS analysis determined that most strains were not related at genetic level. Only 23 (34%) strains were grouped into nine transmission clusters, each one comprising between two and five strains (Fig. 3, Supplementary Fig. S4). These clusters were mostly (91.3%) from Lineage 4. The average genetic distance between the clusters belonging to this lineage was 304.3 nucleotides; while the cluster belonging to Lineage 2 had a considerable minimum distance of 1198 with the nearest cluster of Lineage 4. Likewise, it was evidenced that clusters integrated by strains belonging to the same sublineage showed a lower genetic distance between them compared with the clusters integrated with different sublineages (Fig. 3). All clustered strains came from Lima (n = 21) and Callao (n = 2). Regarding the strains of Lima, 17 belonged to the east zone, 1 to the center (in the area bordering the east zone), 2 to the south zone and 1 to the outskirts (Supplementary Table S4). Regarding the variability of sources of infection, clusters 2 and 6 were composed of strains belonging to a single infection district; while the rest of clusters were integrated by at least two districts. The cluster with more members (cluster 5) was integrated by five strains from the bordering districts of San Juan de Lurigancho and El Agustino. In the same way, cluster 4 was formed by four strains from the geographically close districts of San Juan de Lurigancho, Santa Anita and La Victoria. The remaining clusters were composed of two members, these coming mainly from geographically separated districts (Fig. 3, Supplementary Table S4).

In Peru, the annual number of reported XDR-TB cases is gradually increasing over the years [8]. In our study, we found that most strains belonged to the Euro-American lineage (Linage 4), in agreement with studies that claim that this lineage is the largest in the world and the more prevalent in America continent. In general, the sequenced samples were representatives for XDR-TB affecting the entire country; because, across the years, this form of TB remains highly concentrated in Lima and Callao [35]. Another important finding is that there is significant genomic diversity in our XDR strains and transmission clusters could only be identified in 34% of these.

Concerning the mutations conferring resistance to antituberculosis drugs, currently in Peru, molecular resistance screening for first and second line antituberculosis drugs is performed by the line probe assays (LPAs): GenoType MTBDRplus v2.0 and GenoType MTBDRsl v2.0. However, these assays only concentrate the analysis on genetic hot spots. According to previous studies, the rifampicin resistance is mainly caused by mutations at codons 450, 445 and 435 of rpoB gene ^23–25. We also detected Q432P mutation which has a high confidence grade for rifampicin resistance development and S431R mutation with scarse data according to WHO-NGS Technical Guide, despite having been previously asociaciated with rifampicin resistant phenotypes ²⁶. The rifampicin-discordant strain presented the I491F mutation located outside the RRDR. This mutation was previously reported in Peru ¹⁹ and WHO lists it as a variant with minimum confidence grade ²⁰. The isoniazid resistance was predominantely driven by the S315T mutation of the katG gene followed by c-15t of the inhA promoter. The rare variant inhA g-17t was also present in a low frequency in concordance to previous studies ^15,27. This mutation is indirectly detected by the GenoType MTBDRplus v2.0. Finally, the isoniazid-discordant strain presented the inhA S94A mutation that is known to confer isoniazid resistance in clinical and experimental studies ²⁸. Levofloxacin resistance was mainly driven by mutations occurring at gyrA gene. However, it also was detected in gyrB gene rare mutations E501D and S447F. Some studies establish that this rare mutation only confers resistance to Moxifloxacin or ciprofloxacin ^29–31, although our results indicate that it also confers resistance to levofloxacin and has never been reported before in peruvian strains. Strains that did not harbor mutations in either gyrA or gyrB genes suggests the existence of alternatives mechanisms of resistance like alterations in genes related to efflux pumps as well as DNA mimicry ^32–34. Most mutations occurred in just one gen with exception of a double mutation gyrA G88A and gyrB T500P. this behaviour has been associated with higher MICs and may be associated with a decreased fitness ³⁵, but it could not be determined in this study. All only kanamycin resistant strains exhibited the rrs a1401g mutation which is strongly associated with resistance to high concentrations of kanamycin and amikacin ²⁰. Similarly, the simultaneous resistance to kanamycin and capreomycin is associated with high confidence mutations in the rrs gene (a1401g, c1402t and g1484t) associated with cross resistance to both drugs ²⁰. However, strains with no mutation at rrs, eis or tlyA genes suggest the presence of additional mechanism like alterations in L10 and L12 genes, which have been shown to interact with CAP ³⁶, overexpression of whiB7 and efflux pump genes for KAN ^37,38. This behavior has been previously reported in other studies, including Peru ^39,40.

The analysis of phylogenetic SNPs, together with the spoligotypes, evidenced a predominance of Lineage 4 between the XDR-TB strains analyzed, which is in accordance with studies that claim that this lineage is the largest in the world and the more prevalent in America continent. In addition, the higher proportion sublineages (4.3 and 4.1.2) are considered as worldwide distribution lineages, while Lineage 4.1.1 is considered intermediate distribution worldwide ⁴¹. The high number of strains belonging to Lineage 4 is a characteristic of America, and can be understood due to the colonization of the American continent by European emigrants (founder effect), which is estimated to have occurred approximately between the years 1466 and 1593 ⁴². The low frequency of peruvian XDR-TB strains belonging to Lineage 2 (Beijing family) reveals the recent incorporation of this Lineage into the territory. The same behavior has been reported in countries such as Ecuador and Chile ^43,44. The proportion of strains belonging to the Beijing family (8.8%) is very similar to what was obtained in previous studies conducted in 2012 (9.3%) ⁴⁵ and in 2014 (9.2%) ¹¹. This lineage originated in East Asia and in recent decades it has dispersed globally due to its high degree of association with various virulence factors ⁴⁶. Our study shows that until 2015, the XDR-TB strains belonging to ‘Lineage 2’ could still be geographically restricted to Lima and Callao. Notewhorty, in comparison to spoligotyping, whole genome sequencing allowed to achieve a higher degree of resolution for the differentiation of XDR-TB strains. In this way, several strains with the same spoligotype could be differentiated by genomic analysis to sublineages. Thus, nine strains characterized as belonging to T1 family, were further assigned as belonging to the sublineages: 4, 4.1.1.3 and 4.1.2.1. Moreover, it was evidenced that the only two strains with sublineage 4.1.1.3, present in the geographically remote regions of Arequipa and La Libertad, were separated only by 28 SNPs, suggesting a recent event transmission. Thus, we can appreaciate that the WGS allows to obtain a higher resolution in the classification of MTB strains ⁴⁷.

The small group of strains clustered in recent chains of transmission and the great genetic differences between them suggest that the main mechanism of adquiring a XDR-TB strain is not through by direct contact, but by failures in the individual treatment of less severe forms of tuberculosis. The transmission clusters obtained confirm that XDR-TB strains are concentrated on the districts of Lima that to are more associated with poverty, overcrowding and less access to health systems. These demographic variables have already been well characterized, previously in studies of Peru and the world, but is not being tackled appropriately ⁴⁸.

This study highlights the relevance and utility of performing Whole Genome Sequencing as a high-resolution approach to perform phylogenetics analysis of XDR-TB strains circulating in Peru. We performed Phylogenomic analysis based on both SNPs and spoligotypes evidencing the predominance of Lineage 4 through XDR-TB strains circulating in Peru. Also, transmission analysis indicates that the main mechanism of acquisition of XDR-TB is through failures in the individual treatment of less severe forms of tuberculosis. Finally, the prediction of resistance, determination of transmission groups, and evolutionary analysis can be effectively evaluated using WGS to improve the understanding of XDR-TB dynamics in these settings and provide precise information to improve control measures of TB in Peru.

Sample Collection

MTB strains from hospitals and public health laboratories in Peru are shipped to the Laboratorio de Referencia Nacional de Micobacterias (LRNM), of the Instituto Nacional de Salud (INS), for confirmation and antimicrobial susceptibility testing. A total of 68 XDR-TB strains, according to phenotypic Proportion Method, stored at LRNM were included in this study. All strains correspond to patients with active pulmonary TB recovered between 2011–2015. These strains were randomly selected from the entire country. Approval for the use and processing of preserved strains was obtained from the Institutional Committee for Research Ethics of the INS and the identity of the patients were blinded to the researchers.

Genotypic and Phenotypic Confirmation

All MTB strains were inoculated in Middlebrook 7H9 media (Becton Dickinson GmbH, Heidelberg, Germany) during seven days. Subsequently, 0.2 mL of 7H9 supernatant was transferred to Lowenstein-Jensen (LJ) medium and incubated for a minimum of three weeks to obtain a moderate development. Genotypic confirmations of resistance against rifampicin and isoniazid drugs were performed using the line probe assay (LPA) GenoType® MTBDRplus v2.0 (Hain LifeScience, Nehren, Germany), according to manufacturer's protocol. Phenotypic confirmation was performed using the proportion method (PM) in Middlebrook 7H10 agar BD™ plaques (Becton Dickinson GmbH, Heidelberg, Germany) to asses the resistance against isoniazid (0.2 µg/ml, 1.0 µg/mL), rifampicin (1.0 µg/mL), levofloxacin (1.0 µg/mL), capreomycin (10 µg/mL) and kanamycin (5 µg/mL). All laboratory work related to the culture of live bacteria was performed in the biosafety level 3 (BSL-3) facility of the LRNM.

DNA extraction and Whole Genome Sequencing

Genomic DNA extractions were performed using the GenJET Genomic DNA Purification kit (Thermo Fisher Scientific, Waltham, MA) according to manufacturer's recommendations. Double-stranded DNA concentration was quantified using the Qubit dsDNA HS Assay kit (Thermo Fisher Scientific). Sequencing libraries were prepared using 1 ng of each DNA sample with Nextera XT Library Preparation kit. Whole genome sequencing was done in the INS (Lima, Peru) using Illumina MiSeq platform (Illumina Inc., San Diego, CA, USA) to generate paired-end sequencing reads.

Bioinformatic Analysis

Purity assesment and Quality filtering

Quality evaluation of paired-end reads was performed using fastQC v0.11.9 ⁴⁹. Purity of specific reads for M. tuberculosis complex species was assesed, using a k-mer-based approach, with Kraken2 v2.0.7 ⁵⁰. Verified paired-end reads were filtered with Trimmomatic v0.38 ⁵¹ using default values and a minimum Phred score of 20. Only filtered paired-end reads (95.6% of total raw reads on average) were used for downstream analysis.

Assembly, alignment and Variant Calling

Filtered paired-end reads were mapped against the H37Rv reference genome (GenBank accession number: NC_000962.3) using BWA v0.7.17 ⁵². Duplicate identification and sorting were done with Picard-tools v2.18.25 (http://broadinstitute.github.io/picard). Mapping depth and coverage was determined using samtools v1.9 ⁵³, bedtools v2.29.0 ⁵⁴ and a custom R script. For variant call a local realigment of mapped reads was performed using HaplotypeCaller algorithm, implemented in GATK v3.8 ⁵⁵. A ‘hard filtering’ approach was performed with VCFtools v0.1.16 [23] to select variants with the following criteria: mapping quality ≥ 60, variant depth ≥ 10X and frequency of reads supporting alternate allele ≥ 0.75. Genome positions with missing genotypes and variants identified in repetitive regions (PE, PPE and PE-PGRS families) were also excluded. Finally, selected variants were annotated using SnpEff v4.3T ⁵⁶. Concatenated genome-wide SNPs sequences were generated for subsequent analysis,

Genetic resistance markers

Resistance-associated genes were evaluated to assess phenotypic resistance to rifampicin (rpoB, rpoC, rpoA), isoniazid (katG, inhA, mabA, kasA, furA, ndh, mshA, nat and oxyR-ahpC region), levofloxacin (gyrA, gyrB) and second line injectables (rrs, eis, tlyA) drugs. Variant allelic frequency of at least 0.10 was set for these genes. Genetic variants founded were compared with those reported in the TB Drug Resistance Mutation Database (https://tbdreamdb.ki.se) and confidence was graded based on the technical guide to resistance-associated mutations reported by the WHO ²⁰.

Lineage/sublineage and family determination

MTB lineages and sublineages determination was performed with Kvarq v0.12.2 ⁵⁷ using the set of SNPs proposed by Coll et al ⁵⁸. Families o genotype determination was performed through the in silico detection of 43 unique spacers in the direct repeat (DR) locus using SpoTyping v2.0 [26]. Then, the presence or absence of this spacers were analyzed in the SITVIT2 database (http://www.pasteur-guadeloupe.fr:8081/SITVIT2/) for the determination of the corresponding ‘Spoligo-International-Type’ (SIT).

Evolutionary analysis

A maximum likelihood phylogenomic tree was built from concatenated genome-wide SNPs using RAxML-NG v0.9.0 ⁵⁹, It was used a 25 random and 25 parsimony-based starting trees and 1000 standar non-parametric bootstrap replicates to asses branch support. A general time-reversible (GTR) substitution model was selected based on Akaike’s information criterion using jModelTest2 ⁶⁰. The tree was rooted using a ‘Lineage seven’ strain (SRA ID: ERR181435). Alternative phylogenomic tree was built using additional non-Peruvians 221 XDR-TB strains for a comparison with strains of worldwide distribution (Table S1).

Transmission clusters determination

Genomic transmission clusters were determined using genome-wide SNPs, independently of the epidemiological data. A cutoff value of no more than 10 SNPs distance, pre established for a high prevalence area ⁶¹, was used to group the strains into the same cluster. Maximum SNP distances were obtained from nucleotide pairwise comparisons of all sequenced strains using the Ape package from R v3.6.1 ⁶².

Absence of epidemiological tracing data

ACKNOWLEDGMENTS

We express our thanks to all laboratory staff of the Laboratorio de Referencia Nacional de Micobacterias and to the national network of tuberculosis laboratories, for the routine work in the isolation and identification of different strains that were included in this study. This research was supported by the Instituto Nacional de Salud from Peru and the Programa Nacional de Innovación para la Competitividad y Productividad (INNOVATE-Peru) under the contract 353-PNICP-PIAP-2014.

AUTHOR’S CONTRIBUTIONS

Conceptualization: EDS, LRS, RGG, ANV, ZMP; Formal analysis: EDS, RGG, ZMP; Investigation: EDS, LRS, RGG, ANV, ZMP; Methodology: EDS, RGG, ZMP; Project administration: ZMP; Original draft preparation: EDS, ZMP; Review and editing: EDS, LRS, RGG, ANV, ZMP. The authors read and approved the final manuscript.

COMPETING INTERESTS

The authors declare no competing interests

DATA AVAILABILITY

The data on Whole Genome Sequencing used to support the findings of this study are available from the corresponding author upon request.

World Health Organization. Global Tuberculosis Report 2019. (2019).
Günther, G. et al. Multidrug-Resistant Tuberculosis in Europe, 2010–2011. Emerg Infect Dis21, 409–416 (2015).
Matteelli, A., Roggi, A. & Carvalho, A. C. Extensively drug-resistant tuberculosis: epidemiology and management. Clin Epidemiol6, 111–118 (2014).
Günther, G. et al. Availability, price and affordability of anti-tuberculosis drugs in Europe: a TBNET survey. Eur. Respir. J.45, 1081–1088 (2015).
Mendoza Ticona, A. & Gotuzzo Herencia, E. Tuberculosis extremadamente resistente (TB-XDR): historia y situación actual. Acta Médica Peruana25, 236–246 (2008).
Bonilla Asalde, C. Situación de la tuberculosis en el Perú: current status. Acta Médica Peruana25, 163–170 (2008).
Alarcón, V., Alarcón, E., Figueroa, C. & Mendoza-Ticona, A. Tuberculosis en el Perú: situación epidemiológica, avances y desafíos para su control. Revista Peruana de Medicina Experimental y Salud Publica34, 299–310 (2017).
Rios vidal, J. Situación de Tuberculosis en el Perú y la respuesta del Estado (Plan de Intervención, Plan de intervención).
Barletta, F. et al. Genetic variability of Mycobacterium tuberculosis complex in patients with no known risk factors for MDR-TB in the North-Eastern part of Lima, Peru. BMC Infect. Dis.13, 397 (2013).
Capcha A, L. et al. Perfiles genéticos (IS6110) y patrones de resistencia en aislamientos de M. tuberculosis de pacientes con tuberculosis pulmonar. Lima Sur, Perú. Revista Peruana de Medicina Experimental y Salud Publica22, 4–11 (2005).
Cáceres, O. et al. Characterization of the Genetic Diversity of Extensively-Drug Resistant Mycobacterium tuberculosis Clinical Isolates from Pulmonary Tuberculosis Patients in Peru. PLOS ONE9, e112789 (2014).
Brown, T. S. et al. Genomic epidemiology of Lineage 4 Mycobacterium tuberculosis subpopulations in New York City and New Jersey, 1999–2009. BMC Genomics17, 947 (2016).
Roetzer, A. et al. Whole Genome Sequencing versus Traditional Genotyping for Investigation of a Mycobacterium tuberculosis Outbreak: A Longitudinal Molecular Epidemiological Study. PLoS Med10, (2013).
Walker, T. M. et al. Whole-genome sequencing to delineate Mycobacterium tuberculosis outbreaks: a retrospective observational study. Lancet Infect Dis13, 137–146 (2013).
Walker, T. M. et al. Whole-genome sequencing for prediction of Mycobacterium tuberculosis drug susceptibility and resistance: a retrospective cohort study. Lancet Infect Dis15, 1193–1202 (2015).
Gilchrist, C. A., Turner, S. D., Riley, M. F., Petri, W. A. & Hewlett, E. L. Whole-Genome Sequencing in Outbreak Analysis. Clin Microbiol Rev28, 541–563 (2015).
Niemann, S. et al. Genomic Diversity among Drug Sensitive and Multidrug Resistant Isolates of Mycobacterium tuberculosis with Identical DNA Fingerprints. PLOS ONE4, e7407 (2009).
Schürch, A. C. et al. The tempo and mode of molecular evolution of Mycobacterium tuberculosis at patient-to-patient scale. Infect. Genet. Evol.10, 108–114 (2010).
Solari, L., Santos-Lazaro, D. & Puyen, Z. M. Mutations in Mycobacterium tuberculosis Isolates with Discordant Results for Drug-Susceptibility Testing in Peru. Int J Microbiol2020, (2020).
World Health Organization. Technical guide on next-generation sequencing technologies for the detection of mutations associated with drug resistance in Mycobacterium tuberculosis complex. WHO http://www.who.int/tb/publications/2018/WHO_technical_guide_nextgen_sequencing/en/ (2018).
MINSA. Análisis de la situación epidemiológica en el Perú 2015. (2016).
DGE-MINSA. Análisis de la situación epidemiológica de la tuberculosis en el Perú 2019. (2019).
Agapito, J. et al. Caracterización de las mutaciones en el gen rpoβ asociadas a la rifampicina en pacientes con tuberculosis pulmonar. Revista Peruana de Medicina Experimental y Salud Publica19, 117–123 (2002).
Sandoval, R., Monteghirfo, M., Salazar, O. & Galarza, M. Resistencia cruzada entre isoniacida y etionamida y su alta correlación con la mutación C-15T en aislamientos de Mycobacterium tuberculosis de Perú. Revista Argentina de Microbiología (2019) doi:10.1016/j.ram.2019.03.005.
Farhat, M. R. et al. Rifampicin and rifabutin resistance in 1003 Mycobacterium tuberculosis clinical isolates. J. Antimicrob. Chemother.74, 1477–1483 (2019).
Miotto, P., Cabibbe, A. M., Borroni, E., Degano, M. & Cirillo, D. M. Role of Disputed Mutations in the rpoB Gene in Interpretation of Automated Liquid MGIT Culture Results for Rifampin Susceptibility Testing of Mycobacterium tuberculosis. J Clin Microbiol56, (2018).
Sadri, H., Farahani, A. & Mohajeri, P. Frequency of mutations associated with isoniazid-resistant in clinical Mycobacterium tuberculosis strains by low-cost and density (LCD) DNA microarrays. Annals of Tropical Medicine and Public Health9, 307 (2016).
Shaw, D. J. et al. Disruption of key NADH-binding pocket residues of the Mycobacterium tuberculosis InhA affects DD-CoA binding ability. Scientific Reports7, 1–7 (2017).
Disratthakit, A. et al. Role of gyrB Mutations in Pre-extensively and Extensively Drug-Resistant Tuberculosis in Thai Clinical Isolates. Antimicrob. Agents Chemother.60, 5189–5197 (2016).
Li, Q. et al. Mutation and Transmission Profiles of Second-Line Drug Resistance in Clinical Isolates of Drug-Resistant Mycobacterium tuberculosis From Hebei Province, China. Front Microbiol10, (2019).
Ali, A. et al. Whole Genome Sequencing Based Characterization of Extensively Drug-Resistant Mycobacterium tuberculosis Isolates from Pakistan. PLoS ONE10, e0117771 (2015).
Takiff, H. E. et al. Efflux pump of the proton antiporter family confers low-level fluoroquinolone resistance in Mycobacterium smegmatis. PNAS93, 362–366 (1996).
Louw, G. E. et al. A Balancing Act: Efflux/Influx in Mycobacterial Drug Resistance. Antimicrob Agents Chemother53, 3181–3189 (2009).
Hegde, S. S. et al. A fluoroquinolone resistance protein from Mycobacterium tuberculosis that mimics DNA. Science308, 1480–1483 (2005).
Luo, T. et al. Double mutation in DNA gyrase confers moxifloxacin resistance and decreased fitness of Mycobacterium smegmatis. J Antimicrob Chemother72, 1893–1900 (2017).
Lin, Y. et al. The Antituberculosis Antibiotic Capreomycin Inhibits Protein Synthesis by Disrupting Interaction between Ribosomal Proteins L12 and L10. Antimicrobial Agents and Chemotherapy58, 2038–2044 (2014).
Reeves, A. Z. et al. Aminoglycoside Cross-Resistance in Mycobacterium tuberculosis Due to Mutations in the 5′ Untranslated Region of whiB7. Antimicrob Agents Chemother57, 1857–1865 (2013).
Sowajassatakul, A., Prammananan, T., Chaiprasert, A. & Phunpruch, S. Overexpression of eis without a mutation in promoter region of amikacin- and kanamycin-resistant Mycobacterium tuberculosis clinical strain. Annals of Clinical Microbiology and Antimicrobials17, 33 (2018).
Guio, H., Tarazona, D., Galarza, M., Borda, V. & Curitomay, R. Genome Analysis of 17 Extensively Drug-Resistant Strains Reveals New Potential Mutations for Resistance. Genome Announc2, (2014).
Sowajassatakul, A., Prammananan, T., Chaiprasert, A. & Phunpruch, S. Molecular characterization of amikacin, kanamycin and capreomycin resistance in M/XDR-TB strains isolated in Thailand. BMC Microbiol.14, 165 (2014).
Stucki, D. et al. Mycobacterium tuberculosis Lineage 4 comprises globally distributed and geographically restricted sublineages. Nat Genet48, 1535–1543 (2016).
Brynildsrud, O. B. et al. Global expansion of Mycobacterium tuberculosis lineage 4 shaped by colonial migration and local adaptation. Sci Adv4, (2018).
Jiménez, P. et al. Identification of the Mycobacterium tuberculosis Beijing lineage in Ecuador. Biomédica37, 233–237 (2017).
Meza, P. et al. Presence of Bejing genotype among Mycobacterium tuberculosis strains in two centres of the Region Metropolitana of Chile. Revista chilena de infectología31, 21–27 (2014).
Iwamoto, T. et al. Genetic Diversity and Transmission Characteristics of Beijing Family Strains of Mycobacterium tuberculosis in Peru. PLoS One7, (2012).
Merker, M. et al. Evolutionary history and global spread of the Mycobacterium tuberculosis Beijing lineage. Nat. Genet.47, 242–249 (2015).
Comas, I. & Gardy, J. L. TB Transmission: Closing the Gaps. EBioMedicine34, 4–5 (2018).
Zaman, K. Tuberculosis: A Global Health Problem. J Health Popul Nutr28, 111–113 (2010).
Andrews S. FastQC: a quality control tool for high throughput sequence data. (2010).
Wood, D. E. & Salzberg, S. L. Kraken: ultrafast metagenomic sequence classification using exact alignments. Genome Biology15, R46 (2014).
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics30, 2114–2120 (2014).
Li, H. & Durbin, R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics26, 589–595 (2010).
Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics25, 2078–2079 (2009).
Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics26, 841–842 (2010).
McKenna, A. et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res.20, 1297–1303 (2010).
Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff. Fly (Austin)6, 80–92 (2012).
Steiner, A., Stucki, D., Coscolla, M., Borrell, S. & Gagneux, S. KvarQ: targeted and direct variant calling from fastq reads of bacterial genomes. BMC Genomics15, 881 (2014).
Coll, F. et al. A robust SNP barcode for typing Mycobacterium tuberculosis complex strains. Nat Commun5, 4812 (2014).
Kozlov, A. M., Darriba, D., Flouri, T., Morel, B. & Stamatakis, A. RAxML-NG: a fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference. Bioinformatics35, 4453–4455 (2019).
Darriba, D., Taboada, G. L., Doallo, R. & Posada, D. jModelTest 2: more models, new heuristics and high-performance computing. Nat Methods9, 772 (2012).
Guerra-Assunção, J. et al. Large-scale whole genome sequencing of M. tuberculosis provides insights into transmission in a high prevalence area. eLife4,.
Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics35, 526–528 (2019).

Download PDF

Editorial decision: Major revision
10 Feb, 2021
Reviews received at journal
08 Jan, 2021
Reviewers agreed at journal
07 Jan, 2021
Reviewers invited by journal
07 Jan, 2021
Editor assigned by journal
17 Dec, 2020
Editor invited by journal
10 Dec, 2020
Submission checks completed at journal
10 Dec, 2020
First submitted to journal
04 Dec, 2020

You are reading this latest preprint version

Whole Genome Analysis of Extensively Drug Resistant Mycobacterium Tuberculosis Strains in Peru

Status:

Version 1

Abstract

Figures

Introduction

Results

Patient characteristics

Sequencing and genome assembly

Antimicrobial resistance

* AGC ➡ ACC

Lineages and evolutionary relationships

Transmission cluster determination

Discussion

Conclusions

Methods