The MLVA-5 data for 12,112 S. Typhimurium isolates, which were obtained during routine clinical diagnosis and submitted to the NSW Enteric Reference Laboratory, Pathology West, Sydney for serotyping and MLVA typing in 2007–2014, were analyzed (Sintchenko et al., 2012).

To observe natural VNTR evolution, strains from a novel MLVA type 2-15-8-10-212 (European CDC convention) and its related MLVA types, which emerged from June 2012, were randomly sampled every 2–3 months. Each sampling consisted of strains with MLVA type 2-15-8-10-212 and the MLVA types with one repeat difference which were isolated within 1 month (if available) to represent the diversity of these closely related MLVA types from June 2012 to April 2014. To reduce the impact of geographical diversity, isolates from patients residing only in the eastern suburbs of Sydney were sampled. In total, 28 isolates from MLVA type 2-15-8-10-212 and its related MLVA types were selected for genome sequencing (Table 1, Table S1).

The phenol/chloroform method was used to extract genomic DNA from each strain as described previously (Pang et al., 2013). Genomic DNA was sequenced using the Illumina Genome Analyzer (Illumina) with 300-bp paired-end sequencing. Contigs were de novo assembled using the Velvet version 1.0.8 and VelvetOptimiser (Zerbino and Birney, 2008) and were aligned to the S. Typhimurium LT2 genome (NC_003197, MLVA type 4-13-13-10-211) using progressiveMauve version 2.3.1 (Darling et al., 2010). The coverage of genomes was estimated by Qualimap v2.0 (Garcia-Alcalde et al., 2012). SNP calling was performed as described previously (Pang et al., 2013). A custom script was employed to extract SNPs according to the position on the reference genome. The raw genome sequence data from this study submitted to GenBank is under BioProject No. PRJNA348132.

Maximum likelihood phylogenetic trees were inferred using RAxML 7.2.8 (Stamatakis, 2006). Strain L1880 with MLVA type 2-12-10-10-212 was used as the background strain. Path-O-Gen v1.4 (Rambaut et al., 2016) and BEAST v1.8.2 (Drummond et al., 2012) were used for estimation of genome mutation. The raw estimates in units of substitutions per variable site per year based on a SNP alignment were scaled to genome-wide units of substitutions per site per year by multiplying the estimates by constant k = n/N where n is the number of SNP sites in the alignment and N is the total positions considered for SNP calling (4,487,272; Hawkey et al., 2013). Homoplasy of phylogenetic trees was measured by maximum parsimony analysis in PAUP4.0 (Swofford, 1993). Homoplasy of MLVA was defined as 1 − (K − 1)/M, where K is the minimum number of MLVA mutation events and M is the actual number of VNTR changes along the SNP-tree (Chenal-Francisque et al., 2013).

Three independent replicates of strains L2162 and L2191 representing the two MLVA types, 2-15-8-10-212 and 2-26-6-19-112 respectively, were sub-cultured in Luria-Bertani (LB) liquid culture with two passages per day for a total of 50 passages as previously described (Wuyts et al., 2013). Each replicate was spread onto three LB agar plates and 33 colonies from each replicate were randomly selected. MLVA was performed on these colonies as described previously (Sintchenko et al., 2012).

To evaluate the relationship between SNP and VNTR mutations, pairwise comparisons between isolates were performed by observing SNP difference with a given VNTR locus or repeat difference as done by Eyre et al. (2013). Four pairwise comparisons were performed, including random pairwise comparisons between SNPs and different VNTR locus or repeat difference, pairwise comparisons based on the same time frame within a 1-month window and pairwise comparisons based on the lineages in the SNP tree. The range of SNP differences and the frequency of each SNP difference for all pairwise comparisons of zero, one and two VNTR repeat or locus changes using all VNTR loci were calculated.

Let the mutation rate per genome per generation be denoted by λ_G and the mutation rate of the ith VNTR locus be denoted by μ_i. To characterize the relationship between multilocus VNTR and SNP evolution, we also modeled the time separating two genomes via their most recent common ancestor using the standard coalescent model with an effective population size of N (for an introduction to coalescent theory, see Wakeley, 2008). In coalescent time units where 1 unit is N generations, assuming no homoplasy, the mutation rates generating new SNPs and MLVAs are

Let S and V be the number of SNP and VNTR differences that evolve in these lineages respectively. We made a standard assumption that mutations along a branch occurred according to a homogeneous Poisson process so that over evolution time t the numbers of mutation of each type are distributed as Poisson (λt) and Poisson (μt). Assume that the mutations giving rise to SNPs and VNTRs are independent of each other along a given branch. We computed the joint distribution of the number of SNPs and MLVA differences in a data set as follows. The probability density of the total branch length T in a coalescent for n lineages is

The joint distribution of the number of SNPs and MLVA differences is therefore given by

where p() denotes probability mass. This can be further simplified by changing variables; define x = ((i − 1)/2 + μ + λ)t so that dx/dt = (i − 1)/2 + μ + λ. The integral in the above expression can then be rewritten as

Since υ, s are integer valued, the gamma function Γ(υ + s + 1) can be replaced with the factorial (υ + s)!, so that

The maximum likelihood estimates of λ and μ can be obtained by numerically maximizing the above likelihood equation (Equation 6) with respect to those two parameters. The 95% confidence intervals (CI) for λ and μ were estimated using standard asymptotic theory for likelihood estimators under which standard errors used for the confidence interval are computed using the Fisher Information matrix (Casella and Berger, 2002). Any confidence limit below 0 was set to 0, and each CI includes a Dunn-Sidak correction for testing two hypotheses in a single dataset. The ratio of the two mutation rates (r = λ/μ) was estimated by re-parameterizing the model in terms of r and λ.

We analyzed the allelic distribution of the MLVA data from 12,112 S. Typhimurium isolates from 2007 to 2014. For STTR9 and STTR3, allele 2 and allele 212 were predominant, accounting for 85.5% and 77.3% of the isolates respectively. For STTR9, the next most frequent alleles are 3 and 4 with 9.3% and 4.3% respectively. For STTR3, the next most frequent allele is 211 with 10.0%. There were only a few other alleles at these two loci with low frequencies (Figures 1A,B). STTR3 has a composite of two types of repeats with the unit length of 27 bp and 33 bp, making it more difficult to determine the repeat's composition. Additionally these two loci have been shown to be relatively stable from both in vivo and in vitro passaging experiments (Dimovski et al., 2014). Therefore, they were excluded from further analysis. For the remaining 3 VNTR loci, the distribution of the size of repeat units largely followed a normal distribution, with the median size having the highest frequency by both number of MLVA types and the total number of isolates (Figures 1C–E). We define mid-range repeats as those falling within the second and third quartiles (25–75%) of the repeat distribution. Any repeats outside the mid-ranges are considered to be extreme length repeats. The median number of repeats was 10 and the mid-range was 8–14 repeats for STTR5. Similarly, STTR6 had a median size of 8 and a mid-range of 6–13 repeats while STTR10 had a median size of 11 and a mid-range of 9–12 repeats. We also analyzed the frequency and allele size of the top 10 MLVA types that existed over the 7 years. The top 10 MLVA types vary from year to year from less than 1% to over 20%, which account for 32.6% of total isolates (Table 2). The range of alleles for STTR5, STTR6 and STRTR10 were 7–12, 6–13, and 8 to 13 respectively, all of which were in the mid-range. The predominance of mid-range repeats in the population suggests a direction evolution of the VNTRs toward mid-range repeats.

To confirm the directional mutability of the VNTR loci, an in vitro evolution experiment was carried out using two MLVA types 2-15-8-10-212 and 2-26-6-19-112 represented by L2162 and L2191, respectively. The latter carried an extreme length allele of 26 and 19 repeats for STTR5 and STTR10, respectively. After 50 generations of passaging in three replicates, 34 events were observed in 32 colonies (32/198, 16%), of which 33 mutation events were changes due to single repeats (either insertion or deletion) while one mutation involved change of two repeats (Table 3). The repeat changes fitted the general VNTRs mutation model proposed by Vogler et al. (2006), which predicts that the majority of mutations involve single repeats. Four and three mutation events with one repeat unit deletion were observed in allele 26 and allele 15 (both STTR5), respectively (Figure 2). In contrast, alleles 6 (STTR6), 8 (STTR6) and 10 (STTR10) were all mutated by one repeat unit insertion and observed in 11, seven and five of the 99 colonies for each MLVA type sampled, respectively. These observations suggest that long VNTRs are likely to mutate to shorter ones, and short VNTRs are likely to gradually become longer. Allele 19 (in STTR10) did not exhibit the expected deletion pattern with two insertion events and only one deletion event.

To elucidate the relationship of VNTR change and SNP differences at genome level, we sequenced 28 isolates from 2012 to 2014, from a novel emerging MLVA type, 2-15-8-10-212 and five single locus variant MLVA type at any of the 3 locus, STTR5, STTR6, and STTR10. We chose 2-15-8-10-212 as the main type as the number of this MLVA type (550 isolates) was higher than all of isolates of the other 5 related MLVA types selected and it contained the mid-range size allele for STTR5 (15 repeats) and STTR6 (8 repeats) (Figure S1). We only examined 1 VNTR difference as one allelic change is the most frequent event and more relevant to public health investigations of point source outbreaks (Dimovski et al., 2014; Ahlstrom et al., 2015). We selected 1–6 isolates per MLVA type with matching timeframe as shown in Table S1. All isolates were from eastern Sydney so that the findings are more relevant to local epidemiology as outbreak detection by MLVA is based on spatial and temporal clustering. These isolates were recovered between 20 and 695 days after the initial identification of the first 2-15-8-10-212 isolate (L2162). The VNTR changes and isolation date of the selected 28 S. Typhimurium isolates are shown in Figure 3A.

In total, there were 48 SNPs among the 28 isolates (Table S2). Single nucleotide differences were observed as short as 8 days apart (between isolates L2197 and L2198). Identical isolates were found as far as 151 days apart. Compared with the five strains located near the root of the tree (L2163, L2164, L2166, L2167, and L2197), the isolates showed 1–7 single nucleotide differences. The largest number of single nucleotide differences were found in L2187 with 7 SNPs.

Of the 48 SNPs observed, 43 and 5 SNPs were located in coding regions and integenic regions respectively. The majority of the coding region SNPs (30 out of 43) were non-synonymous SNPs (Table 4). Each of the nsSNPs was located in a different gene, most of the gene functions are related to metabolism. The majority of the nsSNPs were only found in a single isolate, suggesting that they were likely to be transient. nsSNPs present singularly in individual isolates have been observed in other S. Typhimurium isolates with no apparent adaptive values. For example, Octavia et al. (2015b) analyzed the genomic variation of serial S. Typhimurium isolates in human infection related with prolonged carriage and found that the majority of SNPs observed were located in coding regions, most of which were nsSNPs in the genes involved with metabolism. Four nsSNPs were present in 2 or more phylogenetically related isolates including one nsSNP in adrA among 4 isolates (L2175, L2178, L2180, and L2184), one nsSNP in STM1543 among 6 isolates (L2169, L2174, L2175, L2178, L2180, and L2184) and one nsSNP in iroC between 2 isolates (L2179 and L2182). One nsSNP in ccmH was present in all 2-15-8-10-212 isolates except L2171. Three of the four genes are involved in virulence or adaptation. adrA encoding diguanylate cyclase plays a role in biofilm formation (García et al., 2004), STM1543 is involved in swarming motility and virulence (Bogomolnaya et al., 2014), while iroC is involved in iron acquisition (Crouch et al., 2008).

A phylogenetic tree was inferred from the SNPs (Figure 3B). Some isolates were indistinguishable by SNP-based typing but were differentiated by MLVA, and vice versa. For instance, L2178, L2184, and L2180 had the same genome type but different MLVA types. Similarly, two isolates L2163 and L2164 with MLVA type 2-14-8-10-212 and three isolates L2197, L2166, and L2167 with MLVA type 2-15-8-10-212 were found to be genomically identical although the two MLVA types differed only by one VNTR.

The SNPs were mapped on to the branches of the phylogenetic tree (Figure S2). There were no reverse or parallel changes for the SNPs present in more than one isolate, indicating that none of the SNPs was derived from recombination. Maximum parsimony analysis was also performed on the dataset and found that the homoplasy index was zero, further confirming no reverse or parallel changes for the SNPs present in the dataset.

We further examined the relationship of SNPs and VNTRs by pairwise comparisons (Figure 4). The SNP difference per VNTR locus or per VNTR repeat ranged from 0 to 12. The majority of the SNP differences per VNTR repeat fell within 6 SNPs. When we restricted pairwise comparisons to a 1 month window which grouped 28 isolates into 11 monthly intervals (Table S1), the majority of differences were 4 SNPs or fewer. We also examined the SNP differences per repeat along the phylogenetic lineages which have 9 branches (Figure 3B) and compared the isolates at the root with the ones in the respective branches in a pairwise fashion. We found that the range was smaller, ranging from 0 to 7 SNPs. We then compared isolates with 2 VNTR differences. The range of SNP differences was similar to those differing by 0 or 1 (Figure 4), highlighting that there is no relationship between VNTR and SNP differences beyond one VNTR locus difference.

We used the coalescent model to examine the virtual ratio of mutation rates between VNTR changes and genome SNP variation statistically. Using the 28 genomes that contained 48 SNPs and 7 VNTR mutation events (5 mutation events from 2-15-8-10-212 and 2 mutation events from 2-15-9-10-212 in L2184 and L2180), we can obtain a maximum likelihood estimate of per genome per generation mutation rate in effective population size N unit for SNPs and VNTRs of 0.96 (95% confidence interval 0–2.00) and 6.57 (95% confidence interval 1.58–11.56), respectively. The relative rate was estimated to be 1 VNTR change to 6.9 (95% confidence interval 0.65–13.1) SNP changes (Figure S3).

The SNP mutation rate estimated by Path-O-gen (Rambaut et al., 2016) gave a substitution rate of 6.0 × 10⁻⁷ site⁻¹ year⁻¹ or 3 SNPs per genome per year while BEAST (Drummond et al., 2012) gave a substitution rate of 6.83 × 10⁻⁷ site⁻¹ year⁻¹ (3.14–9.62 × 10⁻⁷, 95% HPD) or 3 SNPs genome per year. This rate is similar to a previous estimate from a S. Typhimurium phage type DT135a outbreak in Australia (Hawkey et al., 2013), but faster than that S. Typhimurium outbreaks causing invasive typhoid-like disease in Africa (1–2 SNPs genome⁻¹ year⁻¹; Okoro et al., 2012).

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.