Phylogenomic exploration of the relationships between strains of Mycobacterium avium subspecies paratuberculosis

Abstract

Background: Mycobacterium avium subspecies paratuberculosis (Map) is an infectious enteric pathogen that causes Johne's disease in livestock. Determining genetic diversity is prerequisite to understanding the epidemiology and biology of Map. We performed the first whole genome sequencing (WGS) of 141 global Map isolates that encompass the main molecular strain types currently reported. We investigated the phylogeny of the Map strains, the diversity of the genome and the limitations of commonly used genotyping methods.

Results: Single nucleotide polymorphism (SNP) and phylogenetic analyses confirmed two major lineages concordant with the former Type S and Type C designations. The Type I and Type III strain groups are subtypes of Type S, and Type B strains are a subtype of Type C and not restricted to Bison species. We found that the genome-wide SNPs detected provided greater resolution between isolates than currently employed genotyping methods. Furthermore, the SNP used for IS1311 typing is not informative, as it is likely to have occurred after Type S and C strains diverged and does not assign all strains to the correct lineage. Mycobacterial Interspersed Repetitive Unit-Variable Number Tandem Repeat (MIRU-VNTR) differentiates Type S from Type C but provides limited resolution between isolates within these lineages and the polymorphisms detected do not necessarily accurately reflect the phylogenetic relationships between strains. WGS of passaged strains and coalescent analysis of the collection revealed a very high level of genetic stability, with the substitution rate estimated to be less than 0.5 SNPs per genome per year.

Conclusions: This study clarifies the phylogenetic relationships between the previously described Map strain groups, and highlights the limitations of current genotyping techniques. Map isolates exhibit restricted genetic diversity and a substitution rate consistent with a monomorphic pathogen. WGS provides the ultimate level of resolution for differentiation between strains. However, WGS alone will not be sufficient for tracing and tracking Map infections, yet importantly it can provide a phylogenetic context for affirming epidemiological connections.

Fig. 1

Whole genome SNP-based phylogenetic tree of strains included in this study. a Phylogenetic tree of Map and non-Map strains built using FastTree [52]. Map strains form a single clade and are separated from other mycobacterial species by at least 40,000 SNPs (branches shortened for illustrative purposes). b Maximum likelihood phylogenetic tree of Map strains sequenced as part of this study. The tree was created using RAxML [30], and is based on SNPs identified through mapping to Map-K10 as described in the text. Branches are annotated with the host, country of origin and isolate MAPMRI numbers. Numbers in black to the far right represent the INMV profiles. Previously described lineages are labeled. The dashed box represents strains designated as the Indian bison type. Bootstrap values are shown in Additional file 2: Figure S1

Fig. 2

Type and genomic position of IS1311 element found in the Map sequencing data. The proportion of reads mapping to the IS1311 element with a T or C at position 223 is indicated by the colour gradient of circles displayed on the phylogenetic tree. Asterisks indicate samples where the “TG” deletion was detected [22]. The position of the IS1311 element in the genome was predicted by identifying paired reads where one read maps to the IS1311 element and the other does not, and using the mapped position of the latter. This is displayed horizontally with the genome co-ordinates displayed along the top. The colour intensity of the predicted IS1311 element position represents the relative number of supporting reads detected within the sample, and thus is normalised for differences in depth of coverage between samples

Fig. 3

Investigating the presence of a clock-like signal in the data. The presence of a molecular clock in the dataset was assessed by plotting the root to tip distance of isolates in the phylogeny against isolation date [53] for both the dataset as a whole (a + b) and only Type II isolates (c + d). Both had evidence of a very weak positive signal, as indicated by a low linear regression correlation coefficient (a + c). In order to test the significance of these observations, 99 comparison datasets were produced, in which the isolation dates were permuted on the phylogeny. b and d show histograms of correlation coefficient values from the permuted datasets, with the red line indicating the value from the real data. In both cases the correlation coefficient of the real data is significantly greater than the permuted data at the 0.05 level. Passing this test is a minimal requirement for the application of BEAST analysis (Fig. 4), which assumes the presence of a molecular clock

Fig. 4

Estimates of substitution rate. Coalescent analyses were implemented using the BEAST package (v1.7.5) [31] as described in the text. This was carried out both on the entire data and the Type C isolates alone, using a variety of population models as shown. The mean value represents the mean estimated substitution rate from three independent runs. The confidence intervals represent the maximum and minimum higher posterior densities obtained. As a comparison, the estimated rates are shown alongside a predicted substitution rate of M. tuberculosis in the context of transmission [49]