Longitudinal Outbreak of Multidrug-Resistant Tuberculosis in a Hospital Setting, Serbia

A retrospective population-based molecular epidemiologic study of multidrug-resistant Mycobacterium tuberculosis complex strains in Serbia (2008–2014) revealed an outbreak of TUR genotype strains in a psychiatric hospital starting around 1990. Drug unavailability, poor infection control, and schizophrenia likely fueled acquisition of additional resistance and bacterial fitness–related mutations over 2 decades.

completed following the manufacturer's tissue protocol (from step 4).

Whole-Genome Sequencing (WGS)
DNA library preparation was performed according to the Nextera XT manufacturer instructions (Illumina, http://www.illumina.com). Sequencing was carried out with the Illumina MiSeq and HiSeq system with a minimum average coverage depth of 50-fold. Fastq files have been submitted to the European Nucleotide Archive (ENA);accession numbers are given in Appendix 2. Resulting reads were mapped to the M. tuberculosis H37Rv genome (GenBank accession no. NC_000962.3) using BWA-MEM (2), and mappings refined with the GATK software package (3).
Variants were detected with Samtools (4) and filtered with perl scripts for thresholds of a minimum coverage of 4 reads in both forward and reverse orientation, 4 reads calling the allele with a phred score of >20, and 75% allele frequency. After exclusion of multiple consecutive variant calls (in a 12 bp window), variants in drug resistance associated genes or repetitive regions, the remaining positions that had a clear base call for all strains and matched the above mentioned threshold levels in at least 95% of all strains were considered as valid, and combined in a concatenated sequence alignment. Mutations in drug resistance associated genes were analyzed separately but considering the same thresholds, as mentioned above.

Maximum Likelihood (ML) Phylogenetic Reconstruction and Molecular Clusters
A ML tree was calculated, based on the concatenated sequence alignment using FastTree (5) with a general time reversible (GTR) substitution model and 1,000 resamples. A molecular cluster was defined as >2 strains within a maximum genome wide distance of 5 single nucleotide polymorphisms (SNPs). In case another strain exceeded this threshold but shared a common ancestor with other clustered isolates; we indicated this in Figure 1and in Appendix 2 but did not considered these few cases for the cluster rates.

Comparison of Bacterial Demographic and Molecular Clock Models
For a temporal calibration of phylogenetic trees, we used BEAST 2.4.2 (6) and a concatenated sequence alignment of 6,512 SNPs discriminating all isolates.
Prior to model comparisons within the Bayesian statistical framework, we sought for a proper substitution model using jmodeltest 2.1 (based on maximum-likelihood estimates) (7).
Akaike and Bayesian information criteria (AIC and BIC) were considered for the identification of appropriate substitution models (Appendix 1 Table 2).
To compare a strict versus a relaxed clock model (assuming a constant population size), and different demographic models under the best clock model, we ran each analysis with 300 million steps, sampled 5,000 trees, and discarded 10% as burn-in, resulting in ESS values for all parameters in the hundreds and thousands. Marginal likelihood estimates for each run were  Table 3).

Discrete Phylogeographical Analysis
The spatial and temporal expansion of MTBC strains belonging to the TUR genotype, with a GTR substitution model without invariant sites and no gamma distribution, as well as a tip dating approach (best model for TUR genotype dataset, see below, Appendix 1 Table 3). A Chain length of 300 million ignoring a burn-in of 10%, 5,000 trees sampled, resulted in ESS values for all parameters in the thousands. The temporal and geographic distribution was visualized with SPREAD (11) using Google Earth for location mapping of internal nodes and inferred location changes on the basis of a maximum clade credibility (MCC) tree.

Rationales for Determination of Likely Place of Infection
To identify a possible source/origin for the identified MDR TUR outbreak, i.e., nosocomial transmission versus a possible introduction of a regional-or community-acquired MDR strain type, we compared 2 hypotheses (1). One was of a fast TB progression, assuming that all patients receiving an MDR TB diagnosis at the national center for treatment of all psychiatric patients with concomitant respiratory illnesses Bela Crkva (BC Hospital) have been infected in BC Hospital itself when there was >1 month between admission and diagnosis. If the diagnosis was within the first month following admission at the BC Hospital, the initial hospital or the hometown of the patient was assumed as a likely place of infection (2). The second hypothesis we used was a slow TB progression hypothesis: when MDR TB was diagnosed within the first 2 years after the admission to BC Hospital, we assumed that a patient had been infected elsewhere and was transferred as a latent TB case. The diagnosis after 2 years was considered as a healthcare-associated infection in BC Hospital. The geographic inference of place of infection for both hypotheses was used in the discrete phylogeographic analysis, and marginal likelihoods, inferred by path sampling, were compared as described above.

Statistical Analysis
Categorical data were compared by Pearson's χ 2 test, χ 2 test for trend, and Mann-Whitney U test (or Fisher exact test and Monte Carlo simulation when expected group sizes were less than five). All tests were performed as two-sided tests; p values <0.05 were considered statistically significant. We performed statistical analyses with SPSS version 20 (IBM, https://www.ibm.com/analytics/spss-statistics-software).

Ethics
The study protocol was reviewed and approved by the Ethical Committee, Faculty of Medicine, University of Belgrade, decision number 29/V-14.

DST Results
Phenotypic drug resistance rates, other than those for RIF