Genomic Characterization of Mycobacterium leprae to Explore Transmission Patterns Identifies New Subtype in Bangladesh.

Mycobacterium leprae, the causative agent of leprosy, is an unculturable bacterium with a considerably reduced genome (3.27 Mb) compared to homologues mycobacteria from the same ancestry. In 2001, the genome of M. leprae was first described and subsequently four genotypes (1-4) and 16 subtypes (A-P) were identified providing means to study global transmission patterns for leprosy. In order to understand the role of asymptomatic carriers we investigated M. leprae carriage as well as infection in leprosy patients (n = 60) and healthy household contacts (HHC; n = 250) from Bangladesh using molecular detection of the bacterial element RLEP in nasal swabs (NS) and slit skin smears (SSS). In parallel, to study M. leprae genotype distribution in Bangladesh we explored strain diversity by whole genome sequencing (WGS) and Sanger sequencing. In the studied cohort in Bangladesh, M. leprae DNA was detected in 33.3% of NS and 22.2% of SSS of patients with bacillary index of 0 whilst in HHC 18.0% of NS and 12.3% of SSS were positive. The majority of the M. leprae strains detected in this study belonged to genotype 1D (55%), followed by 1A (31%). Importantly, WGS allowed the identification of a new M. leprae genotype, designated 1B-Bangladesh (14%), which clustered separately between the 1A and 1B strains. Moreover, we established that the genotype previously designated 1C, is not an independent subtype but clusters within the 1D genotype. Intraindividual differences were present between the M. leprae strains obtained including mutations in hypermutated genes, suggesting mixed colonization/infection or in-host evolution. In summary, we observed that M. leprae is present in asymptomatic contacts of leprosy patients fueling the concept that these individuals contribute to the current intensity of transmission. Our data therefore emphasize the importance of sensitive and specific tools allowing post-exposure prophylaxis targeted at M. leprae-infected or -colonized individuals.

Introduction

Mycobacterium leprae and the more recently discovered Mycobacterium lepromatosis (Han et al., 2008) are the causative agents of leprosy in humans as well as animals (Truman et al., 2011; Sharma et al., 2015; Avanzi et al., 2016; Honap et al., 2018; Schilling et al., 2019b; Tió-Coma et al., 2019a). Leprosy is a complex infectious disease often resulting in severe, life-long disabilities and still poses a serious health threat in low- and middle income countries (World Health Organization, 2019). Despite the very limited M. leprae genome variability (Singh and Cole, 2011), the disease presents with characteristically different clinico-pathological forms (Ridley and Jopling, 1966) due to genetically dependent differences in the immune response to the pathogen, resulting in the WHO classification from paucibacillary (PB) to multibacillary (MB) leprosy (Kumar et al., 2017). Notwithstanding the efficacy of multidrug therapy (MDT), approximately 210,000 new cases are still annually diagnosed and this incidence rate has been stable over the last decade (World Health Organization, 2019). Aerosol transmission via respiratory routes is generally assumed to be the most probable way of bacterial dissemination (Bratschi et al., 2015; Araujo et al., 2016). Besides bacterial exposure other risk factors have been shown to be associated with development of leprosy such as genetic polymorphisms (Mira et al., 2004; Zhang et al., 2009; Wang et al., 2015; Sales-Marques et al., 2017), the clinical type of the leprosy index case within a household, immunosuppression (Moet et al., 2004), and nutritional factors (Dwivedi et al., 2019).

Mycobacterium leprae is closely related to Mycobacterium tuberculosis, however, its genome has undergone a reductive evolution resulting in a genome of only 3.27 Mb compared to the 4.41 Mb of M. tuberculosis’ (Cole et al., 2001). Part of the genes lost in M. leprae included vital metabolic activity, causing it to be an obligate intracellular pathogen which cannot be cultured in axenic media that requires support of a host to survive. This poses major limitations to obtain sufficient bacterial DNA for research purposes including whole genome sequencing (WGS). Nevertheless, in 2001 the genome of M. leprae was first published (Cole et al., 2001) leading to the classification of M. leprae into four main genotypes (1–4) (Monot et al., 2005) and subsequently further allocation into 16 subtypes (A–P) (Monot et al., 2009; Truman et al., 2011). The genome of M. leprae contains several repetitive elements such as RLEP which present 37 copies and has been widely applied in molecular diagnostics to specifically detect the presence of this mycobacterium (Donoghue et al., 2001; Truman et al., 2008; Martinez et al., 2011; Braet et al., 2018).

Single-nucleotide polymorphism (SNP) genotyping and WGS are powerful approaches to investigate pathogen transmission as well as bacterial dissemination and evolution through genome characterization (Monot et al., 2005, 2009; Han and Silva, 2014). The limited variation observed in the M. leprae genome permits the reconstruction of historic human migration patterns and the origin of M. leprae (Donoghue, 2019). Over the years, several studies have contributed to the detection and characterization of M. leprae genomes originating from patients all around the world (Monot et al., 2005, 2009; Benjak et al., 2018) as well as from ancient skeletons (Suzuki et al., 2010; Schuenemann et al., 2013, 2018; Mendum et al., 2014; Krause-Kyora et al., 2018), red squirrels (Avanzi et al., 2016; Schilling et al., 2019a; Tió-Coma et al., 2019a), armadillos (Truman et al., 2011; Sharma et al., 2015), non-human primates (Honap et al., 2018), and soil (Lavania et al., 2006, 2008; Turankar et al., 2012, 2014, 2016, 2019; Tió-Coma et al., 2019b). Moreover, skeleton remains have been successfully applied to retrospectively assess whether individuals who contributed to the care of leprosy patients such as the priest Petrus Donders, had developed leprosy (Van Dissel et al., 2019). In the last few years, new tools were developed allowing direct sequencing of M. leprae from various types of clinical isolates (Avanzi et al., 2016; Benjak et al., 2018; Schuenemann et al., 2018). However, these methods were never applied on challenging samples such as slit skin smears (SSS) and nasal swabs (NS) containing a low amount of bacterial DNA compared to skin lesions of patients.

Household contacts of leprosy patients are a high risk group for developing the disease (Richardus et al., 2019), and might serve as asymptomatic carriers contributing to bacterial dissemination. PCR and quantitative PCR (qPCR) are reliable techniques to detect M. leprae DNA and have been proposed as tools for early diagnosis of leprosy, particularly among household contacts of newly diagnosed patients (Gama et al., 2018, 2019). In Brazil, M. leprae DNA has been detected in 15.9–42.4% of healthy household contacts (HHC) in SSS, 9.7–35.2% in blood (Gama et al., 2018; Santos et al., 2018) and 8.9–49.0% in NS (Brito e Cabral et al., 2013; Araujo et al., 2016; Carvalho et al., 2018). Other studies from India, Indonesia and Colombia reported 21% of M. leprae positivity in SSS of HHC (Turankar et al., 2014), 7.8% (van Beers et al., 1994) and 16.0–31.0% in NS (Cardona-Castro et al., 2008; Romero-Montoya et al., 2017).

Detection of host markers, such as serum IgM levels of anti-M. leprae phenolic glycolipid I (PGL-I), represents an alternative approach to diagnose infected individuals (Penna et al., 2016; van Hooij et al., 2017; Barbieri et al., 2019). However, although detection of M. leprae DNA as well as antibodies against PGL-I indicate infection with M. leprae, this does not necessarily result in disease. Thus, these tests alone are not sufficient to identify the complete leprosy spectrum (Spencer and Brennan, 2011; van Hooij et al., 2019).

Bangladesh is a leprosy endemic country reporting up to 3,729 new leprosy cases in 2018 (World Health Organization, 2019). However, M. leprae whole genomes (n = 4) from Bangladesh, have only been described in one study (Monot et al., 2009) in which genotypes 1A, 1C, and 1D were identified. To gain more insight into M. leprae genome variation and transmission routes in endemic areas in Bangladesh as well as the potential role of asymptomatic carriers, we explored the diversity and transmission of M. leprae in four districts of the northwest of Bangladesh. We collected SSS and NS of 31 leprosy patients with a high bacterial load as well as 279 of their household contacts and characterized M. leprae DNA by WGS or Sanger sequencing. The resulting genotypes were correlated to the subjects’ GIS location. Additionally, this is the first study to examine M. leprae DNA detection in comparison to anti-PGL-I IgM levels in plasma measured by up-converting phosphor lateral flow assays (UCP-LFAs).

Materials and Methods

Study Design and Sample Collection

Newly diagnosed leprosy patients (index case, n = 31) with bacteriological index (BI) ≥ 2 and 3–15 household contacts of each index case (n = 279) were recruited between July 2017 and May 2018 (Supplementary Table S1 and Supplementary Data S1) in four districts of Bangladesh (Nilphamari, Rangpur, Panchagar and Thakurgaon). Patients with five or fewer skin lesions and BI 0 were grouped as PB leprosy. Patients with more than five skin lesions were grouped as MB leprosy and BI was determined. The prevalence in the districts where this study was performed was 0.9 per 10,000 and the new case detection rate 1.18 per 10,000 (Rural health program, the leprosy mission Bangladesh, yearly district activity report 2018).

For M. leprae detection and characterization, SSS from 2 to 3 sites of the earlobe and NS (tip wrapped with traditional fiber, CLASSIQSwabs, Copan, Brescia, Italy) were collected and stored in 1 ml 70% ethanol at −20°C until further use. For immunological analysis, plasma was collected (van Hooij et al., 2017, 2018, 2019).

Subjects included in the study were followed up for surveillance of new case occurrence for ≥24 months after sample collection.

DNA Isolation From Slit Skin Smears and Nasal Swabs

DNA was isolated using DNeasy Blood & Tissue Kit (Qiagen, Valencia, CA) as per manufacturer’s instructions with minor modifications. Briefly, tubes containing 1 ml 70% ethanol and SSS were vortexed for 15 s. Slit skin smears were removed and tubes were centrifuged for 15 min at 14000 rpm. Supernatants were removed and buffer ATL (200 μl) and proteinase K (20 μl, Qiagen, provided in the Kit) added. NS were transferred to new microtubes and the microtubes containing the remaining ethanol were centrifuged at 14000 rpm for 15 min. Supernatants were removed and NS were inserted again in the tubes, prior addition of ATL buffer (400 μl) and proteinase K (20 μl). Slit skin smears and NS samples were incubated at 56°C for 1 h at 1100 rpm. Next, AL buffer (200 μl, Qiagen, provided in the Kit) was added and incubated at 70 °C for 10 min at 1400 rpm. Column extraction was performed after absolute ethanol precipitation (200 μl) as per manufacturer’s instructions. To avoid cross contamination tweezers were cleaned first with hydrogen peroxide and then with ethanol between samples.

RLEP PCR and qPCR

RLEP PCR (Donoghue et al., 2001) was performed as previously described (Tió-Coma et al., 2019b). Briefly, the 129 bp RLEP sequence was amplified in 50 μl by addition of 10 μl 5x Gotaq® Flexi buffer (Promega, Madison, WI), 5 μl MgCl2 (25 mM), 2 μl dNTP mix (5 mM), 0.25 μl Gotaq® G2 Flexi DNA Polymerase (5 u/μl), 5 μl (2 μM) forward and reverse primers (Supplementary Table S2) and 5 μl template DNA, water (negative control) or M. leprae DNA (Br4923 or Thai-53 DNA, BEI Resources, Manassas, VA) as positive control. PCR mixes were subjected to 2 min at 95°C followed by 40 cycles of 30 s at 95°C, 30 s at 65°C and 30 s at 72°C and a final extension of 10 min at 72°C. PCR products (15 μl) were used for electrophoresis in a 3.5% agarose gel at 130 V. Amplified DNA was visualized by Midori Green Advance staining (Nippon Genetics Europe, Dueren, Germany) using iBrightTM FL1000 Imaging System (Invitrogen, Carlsbad, CA).

Samples from index cases and a selection of contacts for sequencing were also evaluated by qPCR (Martinez et al., 2009). The mix included 12.5 μl TaqMan Universal Master Mix II (Applied Biosystems, Foster City, CA), 0.5 μl (25 μM) forward and reverse primers (Supplementary Table S2), 0.5 μl (10 μM) TaqMan probe (Supplementary Table S2) and 5 μl template DNA were mixed in a final volume of 25 μl. DNA was amplified using the following profile: 2 min at 50°C and 10 min at 95°C followed by 40 cycles of 15 s at 95°C and 1 min at 60°C with a QuantStudio 6 Flex Real-Time PCR System (Applied Biosystems). Presence of M. leprae DNA was considered if a sample was positive for RLEP qPCR with a cycle threshold (Ct) lower than 37.5 or was positive for RLEP PCR at least in two out of three indecently performed PCRs to avoid false positives.

Library Preparation and Enrichment

A total of 60 DNA extracts were selected for WGS, including 30 from SSS and 30 from NS (Supplementary Figure S1 and Supplementary Data S1). At least one sample from each index case (MB leprosy patient) was selected as well as RLEP positive samples of HHC and MB or PB patients who were household contacts of the index case (selection based on Ct value and household overlap). For 12 subject both SSS and NS samples were selected for WGS. A maximum of 1 μg of DNA in a final volume of 50 μl was mechanically fragmented to 300 bp using the S220 Focused-ultrasonicator (Covaris) following the manufacturer’s recommendations and cleaned-up using a 1.8× ratio of AMPure beads. Up to 1 μg of fragmented DNA was used to prepare indexed libraries using the Kapa Hyperprep kit (Roche) and the Kapa dual-indexed adapter kit as previously described (Benjak et al., 2018) followed by two rounds of amplification. All libraries were quantified using the Qubit fluorimeter (Thermo Fisher Scientific, Waltham, MA), and the fragment size distribution was assessed using a fragment analyzer.

Libraries were target enriched for the M. leprae genome using a custom MYbaits Whole Genome Enrichment kit (ArborBioscence) as previously described (Honap et al., 2018). Briefly, biotinylated RNA baits were prepared using DNA from M. leprae Br4923. A total of 1500 ng of each amplified library was used for enrichment. Each library was pooled prior to enrichment with another library with similar qPCR Ct value. Enrichment was conducted according to the MYbaits protocol with the hybridization being carried out at 65°C for 24 h. After elution, all pools were amplified using the Kapa amplification kit with universal P5 and P7 primers (Roche). All amplification reactions were cleaned up using the AMPure beads (1× ratio).

Illumina Sequencing

Pools were multiplexed on one lane of a NextSeq instrument with a total amount of 20–30 million reads per pool. Some libraries were deep sequenced based on the mapping statistics obtained in the first run.

Raw reads were processed and aligned to M. leprae TN reference genome (GenBank accession number AL450380.1) as previously described using an in-house pipeline (Benjak et al., 2018). A minimum depth coverage of 5 was considered for further phylogenetic analysis.

Sequencing Analysis

Genome comparison was based on analysis of SNPs (analyzed with VarScan v2.3.9; Koboldt et al., 2012) and Indels (analyzed with Platypus v0.8.171; Rimmer et al., 2014) as formerly reported (Benjak et al., 2018). The newly sequenced M. leprae genomes were aligned with 232 genomes available in public databases (Schuenemann et al., 2018; Avanzi et al., 2020). Sites below 80% and above 20% alignment difference were also reported. A comparison to 259 M. leprae genomes (including 27 new genomes) allowed the identification of unique SNPs per index case. Each candidate SNP or Indel was checked manually on Integrative Genomics Viewer (Robinson et al., 2011).

Genotyping and Antimicrobial Resistance by Sanger Sequencing

To further characterize the M. leprae strains for which the whole genome sequence was not obtained, specific primers were designed to perform Sanger sequencing based on unique SNPs (Supplementary Tables S3, S4) of each index case strain. Additionally, Sanger sequencing was performed after amplifying several loci (Supplementary Table S2) to subtype the genomes based on standard the M. leprae classification (Monot et al., 2009; Truman et al., 2011) and to determine antimicrobial resistance to rifampicin (rpoB), dapsone (folP1) or ofloxacin (gyrA). Genotyping by Sanger sequencing was performed to all RLEP PCR positive samples (including samples obtained from leprosy patients and HHC) without a whole genome sequence with a depth coverage of ≥5. PCRs were performed with 5 μl of template DNA using the aforementioned PCR mixes. DNA was denatured for 2 min at 95°C, followed by 45 cycles of 30 s at 95°C, 30 s at 50–58°C and 30 s at 72°C and a final extension cycle of 10 min at 72°C. PCR products were resolved by agarose gel electrophoresis as explained above. PCR products showing a band were purified prior to sequencing using the Wizard SV Gel and PCR Clean-Up System (Promega). Sequencing was performed on the ABI3730xl system (Applied Biosystems) using the BigDye Terminator Cycle Sequencing Kit (Thermo Fisher Scientific). Sequences were analyzed using Bioedit v7.0.5.3.

Anti-PGL-I UCP-LFA

Up-converting phosphor lateral flow assays were performed using the LUMC developed LFA based on luminescent up-converting reporter particles for quantitative detection of anti-M. leprae PGL-I IgM as previously described (van Hooij et al., 2017, 2018, 2019). Plasma samples (n = 308, two samples excluded due to labeling mistake) were thawed and diluted (1:50) in assay buffer. Strips were placed in microtiter plate wells containing 50 μl diluted samples and target specific UCP conjugate (PGL-I, 400 ng). Immunochromatography continued for at least 30 min until dry. Scanning of the LFA strips was performed by LFA strip readers adapted for measurement of the UCP label (UPCON; Labrox, Finland). Results are displayed as the Ratio (R) value between Test and Flow-Control signal based on relative fluorescence units (RFUs) measured at the respective lines. The threshold for positivity for the αPGL-I UCP-LFA was 0.10.

Results

M. leprae Detection in Patients and Healthy Household Contacts

At diagnosis of the index cases and recruitment of contacts in this study out of 279 household contacts 250 presented no signs or symptoms of leprosy or other diseases (HHC), whereas 22 household contacts were diagnosed as PB and seven as MB patients (Supplementary Table S1 and Supplementary Data S1) and therefore were excluded from the HHC group.

Presence of M. leprae DNA was determined by RLEP PCR or qPCR in SSS and NS of leprosy patients and HHC (Figure 1 and Supplementary Data S1): as expected in MB patients with BI 2–6 M. leprae DNA was almost always detectable in both SSS (96.8%) and NS (90.9%). This was much lower in PB and MB patients with BI 0 ranging from 22.2% in SSS to 33.3% in NS. Positivity rates in HHC were not very different from those observed for PB and MB patients with BI 0, with 12.3% positive samples in SSS and 18.0% in NS. Showing a similar M. leprae carriage between HHC and patients with BI 0. Moreover, the overall Ct range was lower for SSS [16.3–37.1] compared to NS [20.1–39.4] showing that SSS contained more M. leprae DNA and is a preferred sample for its detection (Supplementary Data S1).

FIGURE 1

Study design, RLEP positivity and genotyped samples. Flow diagram providing an overview of the subjects recruited for this study. Slit skin smears (SSS) and nasal swabs (NS) collected per group; healthy household contacts (HHC), paucibacillary (PB) or multibacillary (MB) patients with BI 0, and MB patients with a bacteriological index (BI) 2–6. MB patients with BI 1 were not diagnosed within the course of this study. DNA was isolated from SSS and NS and screened for M. leprae DNA by RLEP PCR. Samples were genotyped by Sanger sequencing (Monot et al., 2009; Truman et al., 2011) or Whole Genome Sequencing (Benjak et al., 2018). Percentages of the samples positive for RLEP PCR and genotyped are shown.

HHC (n = 250) were followed up clinically for ≥24 months after sample collection and four of them developed leprosy within the first year. RLEP PCR performed on DNA isolated before disease occurrence showed a positive result from SSS in one patient (5 months before diagnosis) and a positive result from NS in another (8 months before diagnosis). All of the new cases developed PB leprosy with BI of 0 and three were genetically related to the index case (parent and child of index case H03 and second degree relative of index case H30) and one was the spouse (index case H10).

Genome Typing and Antimicrobial Resistance

Mycobacterium leprae genomes of SSS and NS were genotyped by WGS or Sanger sequencing. A total of 60 samples (30 SSS and 30 NS) from MB and PB leprosy patients as well as HHC were selected for WGS with an RLEP qPCR Ct ranging from 16.2 to 37.2 (Supplementary Data S1). A total of 27 samples from 21 subjects (21 SSS and 6 NS) passed the library quality check and were successfully sequenced with a coverage ≥5 (Supplementary Figure S1 and Supplementary Table S5). The limiting Ct value was 26.2 for SSS and 24.2 for NS.

On applying the genotyping system described by Monot et al. (2009), Truman et al. (2011), the following genotypes were found for these 21 subjects: 1A (n = 5), 1B (n = 4), 1C (n = 3), and 1D (n = 9). Interestingly, the four newly sequenced 1B genotype strains do not cluster with the two previously described 1B strains from Yemen and Martinique (Figure 2). Instead, they form a new cluster in the phylogenetic tree located between genotypes the 1A and 1B, which we refer to as 1B-Bangladesh (Figure 2, blue and Supplementary Data S1). Using Sanger sequencing, the M. leprae strain for eight additional individuals were determined as 1A (n = 4) or 1D (n = 4). Three subjects, including two NS samples from HHC, carried genotype 1 but subtype could not be established due to lack of amplification of the subtyping loci (Supplementary Data S1).

FIGURE 2

Phylogeography of M. leprae strains. Maximum parsimony tree of 259 genomes of M. leprae built in MEGA 7. Support values were obtained by bootstrapping 500 replicates. Branch lengths are proportional to nucleotide substitutions. The tree is rooted using M. lepromatosis. The strains from Bangladesh are shown in red and their exact organization in the tree is shown in the two zoomed sections of the genotypes 1A-B and 1D. Strains with an A at SNP61425 in the esxA gene are shown in green. The specific 1B-Bangladesh genotype/cluster of Bangladesh strains is shown in blue.

The SNP used to differentiate genotype 1C (A61425G; Met90Thr, mutated in genotypes 1D and 2–4) is located at esxA. In contrast to previous observations (Monot et al., 2009; Truman et al., 2011), we found that this position is not phylogenetically informative as it is also found unmutated (A; Met) in strains from the genotype 3I and 2E (Figure 2, green and Supplementary Data S2). Moreover, the 1C strains clustered in the middle of the 1D group suggesting that the previously described genotype 1C is part of the 1D genotype.

Finally, antimicrobial resistance was assessed in all genotyped strains either by WGS or Sanger sequencing. The latter was successful on 18 samples for rpoB, five sample for folP1 and 15 samples for gyrA (Supplementary Data S1). None of the strains with a complete genome harbored drug-resistance mutations. One NS sample containing a missense mutation in the rpoB gene (Ser456Thr) in 50% of the sequences potentially leading to antimicrobial resistance (World Health Organization, 2017) was identified by Sanger sequencing. Moreover, although not causing resistance, up to two silent mutations in three different positions of the rpoB gene relevant for antimicrobial resistance (432, 441, and 456) were also observed in several subjects.

Distribution and Possible Transmission of M. leprae Genotypes

The most prevalent M. leprae genotype in the studied area of Bangladesh is 1D, found in 55% of the individuals (n = 16, Table 1 and Supplementary Data S1), followed by 1A in 31% (n = 9), and 1B-Bangladesh in 14% (n = 4). Genotype 1D is the most widely distributed throughout the whole area studied (Figure 3, blue and purple), whilst genotypes 1A and the here identified genotype 1B-Bangladesh are only observed in the eastern area (green and orange, respectively). The latter genotype was found in four individuals: two from the same household and two unrelated subjects residing 56, 51, and 11 km from each other. However, due to privacy regulations on patient information to third parties it could not be established whether subjects in different households had had contact with any of the others.

TABLE 1

M. leprae genotypes identified in Bangladesh.

Genotype	Number of individual	%
1A	9	31.0
1B-Bangladesh	4	13.8
1D	13	44.8
1D-esxA	3	10.4
1*	3

FIGURE 3

Distribution of M. leprae genotypes in Bangladesh. Map of Bangladesh including markers indicating the residence of every subject with at least one sample genotyped for M. leprae (A), and zoomed into the area of interest (B). Each marker indicates an individual for whom M. leprae genotype was determined, either from slit skin smear, nasal swab or both samples. Genotype 1A is shown in green, 1B-Bangladesh in orange, 1D in blue, 1D-esxA in purple and 1* in white. 1D-esxA is 1D subtype containing an A at SNP61425 in the esxA gene, formerly grouped as 1C (Monot et al., 2009; Truman et al., 2011). 1* are samples with genotype 1 for which the subtype could not be determined. The figure was drawn in R (v3.4.3) with the package leaflet (v2.0.2) using maps from Esri – National Geographic with permission. Scale Not Given. “National Geographic World Map.” December 13, 2011. http://www.arcgis.com/home/item.html?id=b9b1b422198944fbbd5250b3241691b6 (September 2, 2019).

In a total of four households the same M. leprae genotype was detected in two individuals (Supplementary Data S1). In the first household, both subjects were MB patients and WGS showed no genetic variation between both patients’ genomes (RB001 and RB003, 1B-Bangladesh genotype, Supplementary Data S2). In the second household with two MB patients, the M. leprae whole genome was only obtained from the index case but the same genotype, 1A, and a strain-specific SNP of the index case (Supplementary Tables S3, S4) was also identified by Sanger sequencing in the other patient (RB182 and RB266). In the last two households, the genotype of strains from both MB index cases’ were determined by WGS (RB030, genotype 1D) and, by Sanger sequencing (RB065, genotype 1D-esxA), while the M. leprae genotype 1 was located in the NS of both HHC but no further subtyping was possible.

Comparison of M. leprae Genomes From SSS and NS

Mycobacterium leprae whole genomes of six patients were successfully recovered from both SSS and NS. The M. leprae genotypes obtained in each subject were in agreement between the two samples (Table 2). Genomic comparison showed no differences between DNA from SSS and NS for two patients: RB001-RN001 (genotype 1B-Bangladesh) and RB048-RN059 (genotype 1D-esxA, Figure 2 and Supplementary Data S2).

TABLE 2

Intraindividual M. leprae genomic differences.

Samples	Genotype	Mutation and amino-acid change	Genomic region	NS		SSS		Other genomes with similar position mutated (variant,% reads, genotype)
				% reads	Aligned reads	% reads	Aligned reads
RB073-RN084	1A	T1824441C;Gly56Asp	ml1512	17%	24	-	-	-
RB053-RN022	1D	G1823127A; Ser494Leu	ml1512	35%	26	-	-	Ml10-98 (Ser494Ala, 95%, 4N)
		1823613_1823614insC;Asp332fs	ml1512	20%	21	-	-	ARLP-23 (Asp332fs, 80%, 2E)
RB074-RN095	1B- Bangladesh	G1823098A; Leu504Phe	ml1512	-	-	92%	92	-
		G660474C; Val252Leu	metK	75%	12	100%	92	-
		C2116695A; Pro100Thr	ml1750	48%	21	-	-	SM1, (Pro100Ser, 100%, 4N/O)
		A2116670G; Gln108Arg	ml1750	20%	20	-	-	LRC-1A (Gln108His, 100%,1A) Md0536 (Gln108Arg 90%, 1D-Mada)
		C2116490T; Arg168His	ml1750	19%	16	-	-	Br14-3 (5 Arg168Cys, 1%,3I) Arg168His: Br2016-17 (22%, 3I); Comore-3 (36%; 1D-Mada); Ml9-81 (29%; 4N)
RB069-RN165	1B- Bangladesh	C95231T; Leu34Phe	glpQ	29%	76	-	-	-
		C2121552T; Val226Ile	ml1752	15%	94	-	-	-

In a third patient (RB073-RN084, genotype 1A), both strains were identical except that in the NS strain 17% of 24 reads in ml1512 harbored a T1824441C (Gly56Asp) (Table 2). Interestingly, ml1512 which encodes a ribonuclease J is one of the most mutated genes among all M. leprae strains (Benjak et al., 2018) and mutations at this gene were also observed in two different patients: in the NS of RN022-RB053 (genotype 1D) 35% of 26 reads had a mutated allele (G1823127A; Ser494Leu) and 20% of 21 reads had an insertion of a C at position 1823613 probably leading to a deleterious frameshift; in the SSS of RB074-RN095 (genotype 1B-Bangladesh) 92% of 92 reads presented a missense mutation (G1823098A; Leu504Phe). Interestingly, RB074 harbored a G660474C mutation in metK, a probable methionine adenosyl-transferase, which was also found in 75% of 12 reads of the NS and is uniquely found in this subject’s M. leprae genomes. Additionally, RN095 also displayed mutations at several positions in ml1750 (a putative nucleotide cyclase): 48% of 21 reads had C2116695A (Pro100Thr), 20% of 20 reads had A2116670G (Gln108Arg) and 19% of 16 reads had C2116490T mutation (Arg168His). These positions were partially or totally mutated in other strains from different genotypes: SM1 (100% Pro100Ser; genotype 4), Ml9-81 (Mali, 30% Arg168His; genotype 4N) and Md05036 (Madagascar, 90% Gln108Arg, genotype 1D-Mada) (Benjak et al., 2018; Avanzi et al., 2020).

The patient with the M. leprae strains that were the most genetically different between the NS and SSS carried the genotype 1B-Bangladesh (RB069 and RN165). The NS strain had a mixed population in glpQ (29% of 76 reads C9231T, Leu34Phe) and ml1752 (15% of 94 reads C2121552T, Val226Ile). These genes encode a glycerophosphoryl diester phosphodiesterase, a putative nucleotide cyclase and a conserved hypothetical protein. Notably, ml1752 is also one of the most hypermutated genes in M. leprae (Benjak et al., 2018).

For 11 patients a whole genome sequence was recovered only from SSS but Sanger sequencing was successfully performed to identify the subtype in NS. The same subtype observed in SSS was also found in the NS of these 11 patients. Moreover, unique M. leprae SNPs identified in the genomes of the SSS (Supplementary Tables S3, S4) were also detected in seven of the genomes of the NS of these patients (Supplementary Data S1).

Combining Host and Pathogen Detection

Anti-PGL-I IgM levels were determined in plasma of 308 subjects. All MB patients with BI 2-6 (n = 33) showed high levels for anti-PGL-I IgM (Table 3) in line with the general consensus (Geluk et al., 2011; van Hooij et al., 2017). Out of the patients (both MB and PB) with BI 0 (n = 27), nine (33.3%) were positive for anti-PGL-I IgM. Similarly, 36.8% of HHC showed positivity (n = 92). From these 92 positive individuals, 70 were neither positive for SSS nor NS RLEP PCR (Supplementary Data S1).

TABLE 3

Anti-PGL-I IgM positivity.

Genotype	Number of positive individual	% of positivity
MB patients BI 2-6 (n = 33)	33	100.0
Patients BI 0 (n = 27)	9	33.3
Healthy household contacts (n = 250)	9	36.8

Of the four contacts who developed leprosy within the first year after sample collection, two were positive for anti-PGL-I IgM whilst negative for RLEP PCRs 10 and 12 months before diagnosis. Since the two other subjects had a positive RLEP PCR in SSS or NS 5 or 8 months before diagnosis, it can be concluded that all of the new cases showed positivity either for host- or pathogen-associated diagnostics 5-12 months before developing disease.

Individual anti-PGL-I levels were compared to RLEP Ct values in SSS and NS samples (Figure 4), showing an expected negative correlation between anti-PGL-I ratio and Ct value since both values are associated with BI. A subtle difference can be observed in the correlation between anti-PGL-I IgM levels and RLEP Ct if the qPCR was performed on either SSS or NS DNA, with a coefficient of determination (R2) 0.73 and 0.69, respectively.

FIGURE 4

Correlation of IgM antibodies against PGL-I to Ct of RLEP qPCR. Quantified levels of pathogen DNA (qPCR) and host immunity were correlated for samples selected for qPCR analysis based on RLEP positivity in multiple individuals in one household. Each dot represents a sample from one individual; leprosy patients are indicated in black, and healthy household contact in blue. Anti-PGL-I antibody levels were measured by up-converting reporter particles lateral flow assay specific for M. leprae PGL-I IgM antibodies (αPGL-I UCP-LFA) using the Ratio (R) of the Test (T) and flow control (FC) lines as units. Ratios of ≥0.10 were considered positive as indicated by the red dashed line. RLEP cycle threshold (Ct) values are indicated on the x-axis and were measured by qPCR to detect M. leprae DNA in slit skin smears (SSS, left) and nasal swabs (NS, right). Undetermined Cts are depicted as Ct 40.

Discussion

In this study we investigated M. leprae transmission patterns in Bangladesh by detecting and sequencing M. leprae DNA derived from SSS and NS of patients and their household members. Our data represents the first report of M. leprae DNA detection in HHC from Bangladesh. We observed moderate positivity in HHC which was similar to positivity of leprosy patients with BI 0. A new genotype, 1B-Bangladesh, was sequenced and we showed that the previously described 1C genotype is part of the 1D group. Additionally, a negative correlation between RLEP Ct values indicating the amount of M. leprae DNA and anti-PGL-I IgM levels was observed.

Mycobacterium leprae DNA detection frequency in HHC from Bangladesh (12.3% in SSS and 18.0% in NS) was in line with previous studies conducted in several hyperendemic areas of Brazil, Colombia and Indonesia (van Beers et al., 1994; Cardona-Castro et al., 2008; Brito e Cabral et al., 2013; Romero-Montoya et al., 2017; Carvalho et al., 2018; Gama et al., 2018). In India higher positivity (21.0%) in SSS of HHC was reported (Turankar et al., 2014) whereas in two Brazilian studies from Uberlandia, up to 42.4% positivity in SSS (Santos et al., 2018) and 49.0% in NS (Araujo et al., 2016) were observed. Three factors may limit the translation of these high positive results from India and Brazil to our study: (i) the sample sizes of the Indian (Turankar et al., 2014) and one of the Brazilian studies (Araujo et al., 2016) were smaller (n = 28 and n = 104, respectively vs. n = 250 HHC in this study); (ii) we conducted a more stringent approach by testing the samples in three independent PCRs; and (iii) the epidemiology and incidence of MB cases in India and Brazil differ from the studied area in Bangladesh where MB leprosy cases occur less frequently than PB and also usually display a low BI (Richardus et al., 2017; van Hooij et al., 2019).

Mycobacterium leprae DNA in the nose does not indicate disease but (transient) colonization whilst presence of M. leprae in SSS indicates infection. Thus, the higher RLEP PCR positivity in NS compared to SSS in patients with BI 0 and HHC likely represents the (virtual) absence of bacteria causing infection in these individuals despite colonization.

A longitudinal study conducted in Brazil (Manta et al., 2019), investigated SSS from 995 HHC by qPCR including follow-up for at least 3 years with occurrence of five new cases. The authors reported 20% qPCR positivity in HHC representing future new cases compared to 9% in HHC without disease. However, this difference was not significant. In line with that study, we found that M. leprae DNA detection was slightly higher (25% vs. 18% in NS and 25% vs. 12% in SSS) in contacts who developed disease compared to those who did not. Additionally, we determined anti-PGL-I IgM levels, which correlated well with Ct qPCR values. Notwithstanding this correlation, serology provided added value: when positivity in any of the three techniques was considered (NS PCR, SSS PCR, or anti PGL-I), all of the contacts (n = 4) who developed leprosy within the first year after sample collection, were identified. In agreement with this, a combination of host and pathogen markers was previously integrated in a machine learning model using qPCR and serological data (antibodies against LID-1 or ND-O-LID) (Gama et al., 2019) to identify prospective leprosy patients among contacts leading to an increased sensitivity in diagnosis, particularly in PB leprosy. It is of note that in our study, three of the four contacts who developed leprosy were genetically related to the index cases in their households, stressing the previously described role of genetic inheritance in the development of leprosy (Mira et al., 2004; Zhang et al., 2009; Wang et al., 2015; Sales-Marques et al., 2017; Uaska Sartori et al., 2020). For this reason, the association between leprosy and the genetics of this Bangladeshi population is currently being studied.

Genotype 1 was identified in all the M. leprae genomes retrieved from Bangladesh, consistent with previous data from Monot et al. (2009). In Bangladesh, leprosy was likely introduced through the southern Asian route (genotype 1) leading to the spread of M. leprae into the Indian subcontinent, Indonesia and the Philippines (Monot et al., 2009; Benjak et al., 2018). Subtype 1D was predominantly present in Bangladesh but in addition we detected 1A and identified a new 1B-Bangladesh genotype. The presence of multiple subtypes of M. leprae genotype 1 in Bangladesh is in line with previous studies in South Asian countries such as India, Nepal, Thailand, Indonesia and Pakistan (Monot et al., 2009; Benjak et al., 2018). The new 1B-Bangladesh genotype is thus far restricted to Bangladesh and two of the four individuals carrying this strain were part of the same household whilst the other two did not have any relationship with each other and were located in different areas with a distance of up to 56 km between them. This suggests that this genotype could be a common subtype in Bangladesh although additional studies are required to confirm this. Thus, it is of interest to include the 1B-Bangladesh SNP specific primers in future epidemiological studies, particularly in other (neighboring) Asian countries such as India where genotype 1 is widely established (Monot et al., 2009).

In contrast to the general belief (Monot et al., 2009; Truman et al., 2011), we observed that subtype 1C does not form an independent subtype but actually belongs to subtype 1D. SNP61425 used to distinguish genotypes 1A–C is located at esxA encoding the virulence factor ESAT-6 (Monot et al., 2009). The Esx protein family also revealed high diversity in the more pathogenic mycobacterium, M. tuberculosis (Uplekar et al., 2011), and is involved in host-pathogen interaction. Of note is that ESAT-6 (ML0049) is a potent T-cell antigen (Geluk et al., 2002, 2004), thus mutations in esxA gene might indicate drift due to immune pressure potentially explaining the occurrence of mutations at SNP61425 in different genotypes.

In a recent survey in 19 countries during 2009–2015 (Cambau et al., 2018), 8% of the cases presented mutations resulting in antimicrobial resistance and resistance to up to two different drugs was detected. In our study, which is the first investigating M. leprae drug resistance in Bangladesh, we detected no resistance by WGS, however, a partial missense mutation in the codon for Ser456 of the rpoB gene potentially leading to rifampicin resistance (n = 1) was observed by Sanger sequencing. This could be the result of a mixed infection or an emerging mutation of the M. leprae strain occurring in the patient. Silent mutations in the rpoB gene were detected in several locations, which indicates that mutations do occur, and this may eventually lead to missense mutations conferring antimicrobial resistance. However, drug resistance is not only induced by genetic mutations in drug targets, efflux systems resulting in antimicrobial resistance have also been described for M. leprae (Machado et al., 2018). This mechanism of drug resistance is unnoticed in genomic tests and needs to be further investigated for leprosy especially in the light of the huge efforts recently initiated and WHO-endorsed for post-exposure prophylaxis (PEP) using antibiotic regimens (Barth-Jaeggi et al., 2016; Mieras et al., 2018; Richardus et al., 2019).

Despite our finding that NS samples were more frequently positive for M. leprae DNA, recovery of M. leprae whole genomes from SSS has proven to be more successful than from NS. This is due to the higher number of bacteria in SSS of patients. However, the importance of genotyping NS as well as skin biopsies or SSS to better understand transmission has been previously discussed (Fontes et al., 2017), as the nasal respiratory route remains one of the most plausible modes of infection (Bratschi et al., 2015; Araujo et al., 2016). In a recent study, skin biopsies and NS of patients were compared by VNTR typing and the authors found that out of 38 patients, differences between SSS and NS in seven loci were observed in 33 patients (Lima et al., 2016). Although the M. leprae genomes from SSS and NS analyzed in our study were almost identical, we observed that genomes obtained from NS harbored more mutations, especially in previously reported (Benjak et al., 2018) hypermutated genes. This could be an indication of in-host evolution in the nasal mucosa, mixed infection or mixed colonization. Thus, it may imply that colonization occurred with two different strains causing a co-infection or that one is present, likely from a later colonization, but does not cause the disease.

The presence of mixed infections emphasizes once more the importance of monitoring asymptomatic carriers, who may contribute to the spread of the pathogen. Therefore, providing PEP only to the (close) contacts of leprosy patients might not be sufficient to stop transmission. Instead, an approach including the entire community but targeting only individuals testing positive for M. leprae DNA or host immune markers associated to M. leprae infection, would represent a preferred strategy for PEP.

Author’s Note

This manuscript has been released as a pre-print at medRxiv (Tió-Coma et al., 2020).

Data Availability Statement

Sequence data are available from the NCBI Sequence Read Archive (SRA) under the bioprojects PRJNA605605 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA605605) and PRJNA592722 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA592722), biosamples SAMN14072760-775 and SAMN13438761-771.

Ethics Statement

The studies involving human participants were reviewed and approved by the National Research Ethics Committee (BMRC/NREC/2016-2019/214). Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

Author Contributions

AG, MT-C, and JR conceived the study. MT-C and JR contributed to the data curation. MT-C, CA, AH, and AB contributed to the formal analysis. AG and JR acquired funding. MT-C, CA, EV, LP, and AH investigated the study. MK, KA, PC, and SC contributed to the resources. AG supervised the study. MT-C wrote the original draft. MT-C and AG wrote, reviewed, and edited the manuscript. All authors reviewed, discussed, and agreed with the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.