Using whole genome sequencing to investigate transmission in a multi-host system: bovine tuberculosis in New Zealand

Sampling and isolate preparation

As part of the routine bTB surveillance in New Zealand, any cattle or wildlife suspected of M. bovis infection undergo a post-mortem examination, and if lesions are discovered a selection are investigated using culture and strain typing. Conventional tests (described in de Lisle et al. [24]) were used to positively identify M. bovis infection. Isolates were REA typed according to previously described methods [16, 25] and cultures were frozen and stored in the strain archive at AgResearch Ltd. Isolates from the archive were selected to provide a representative sample of the M. bovis population circulating in cattle and wildlife across New Zealand between 1985 and 2013 (Fig. 1). To create this representative sample, groups of isolates, from cattle and wildlife, of the same or closely related REA type from the same geographical region were selected from the Central North Island region, and the West Coast and Northeast regions of the South Island. The groups were selected to include all of the most frequently isolated REA types.

Selected isolates were re-cultured at AgResearch Ltd. to generate DNA (Deoxyribonucleic Acid) extracts for WGS. Frozen culture stocks were grown to mid log phase (OD600 = 0.4–0.8) in 5 ml of Tween/albumin broth [26] and sub-cultured into 100 ml of the same media for 5 to 11 weeks until the cultures reached stationary phase. Cell cultures were heat killed and stored at −20 °C. Bacterial DNA was specifically separated from the other cellular components with a high salt hexadecyl trimethyl ammonium bromide (CTAB) extraction. DNA was then extracted with a chloroform/isoamyl alcohol method and quantified using Invitrogen Qubit fluorometry.

A selection of 90 DNA isolates were sequenced at the University of Glasgow Polyomics Facility using an Illumina MiSeq platform that produced 2 × 300 bp paired end reads per isolate. 17 additional isolates were sequenced at New Zealand Genomics Ltd. on a MiSeq platform that produced 2 × 250 bp paired end reads. The remaining isolates (n = 204) were sequenced at the Wellcome Trust Sanger Institute using an Illumina HiSeq that produced 2 × 100 bp paired end reads.

Processing sequencing data

The raw reads for each isolate were examined using FASTQC (v0.11.2 - Andrews [27]) to identify poor quality ends that were then trimmed using PRINSEQ (v0.20.4 - Schmieder & Edwards [28]). If adaptor sequences were present these were removed using TRIMGALORE (v.0.4.1 - Krueger [29]). Each isolate’s trimmed reads were aligned to the M. bovis reference genome, AF2122/97 [30], using the freely available Burrows-Wheeler Alignment tool [31, 32]. The mean coverage (sites with Read Depth (DP) ≥ 20) for the isolates was 99% (2.5% Lower: 96.9, 97.5% Upper: 99.8).

Site information across the isolates was collated to allow the quality of individual sites to be assessed. Sites that fell within Proline-Glutamate (PE) and Proline-Proline-Glutamate (PPE) genes or annotated repeat regions were removed (Sampson [33]). Thereafter only Variant Positions (VPs), sites for which at least one of the isolates showed variation against the reference genome, were retained.

High quality sites were selected for subsequent analyses based on the Mapping Quality (MQ), High Quality base Depth (HQDP) and Read Depth (DP). Filters were designed using MQ, HQDP, DP, in addition to the support for the allele called (SUP), the site coverage across the isolates (COV), and the number of positions that separated SNPs (PROX). A filter sensitivity analysis was conducted in order to establish an optimal filter combination (See Additional file 1: Supplementary text 1.1). The following filters were selected: MQ ≥ 30, HQDP ≥ 4, DP ≥ 30, SUP ≥ 0.95, COV ≥ 0.7, and PROX = 10.

Isolate selection

An early examination of the WGS data revealed that, although most isolates with the same Restriction Endonuclease Analysis (REA) type were very similar, several isolates were quite distinct from the others with the same REA type. These “outliers” were further investigated to determine whether they were mislabelled. Although it was not possible to re-examine these isolates with the REA typing method, potentially mislabelled samples were further examined with Variable Number Tandem Repeat (VNTR) assays. Specific REA types are known to be associated with specific VNTR types. VNTR assays were conducted (described in [15]) using a subset of VNTR loci that were likely to discriminate the isolates in question. For controls, a selection of isolates with similar sample numbers to the questionable isolates was also re-examined. If the determined VNTR types differed from what would have been expected based on previous analyses of these types [15], the isolate was considered to have been mislabelled. 15 of the 28 (14 suspects and 14 controls) isolates examined had VNTR loci that differed from what was expected. These 15 mislabelled isolates were removed from any further analyses (See Additional file 1: Supplementary text 1.2), leaving 296 isolates for further investigation.

Using the 296 isolates, a maximum likelihood phylogenetic tree was constructed in the program PHYLIP (v3.695 - Felsenstein [34]) and rooted using the program PATH-O-GEN (v1.4 - Rambaut [35]). For each isolate the sampling location (including latitude and longitude) and year (of sample submission), sampled species, and REA type were available. Using the maximum likelihood tree and the available sampling information, a selection of spatially and temporally associated isolates were chosen from within clade 1 (Fig. 1a).

Although a large number of isolates were available for the current analyses, these isolates fell within highly distinct clades. Isolates from a single clade were selected to ensure a relatively recent common ancestor to the isolates analysed, and limit the effects of biases introduced by examining genetically distinct groups. Unique pairs of cattle and wildlife isolates were chosen from within strict spatial (40 km) and temporal (+/− 3 years) limits to reduce the impact of potential temporal and spatial sampling biases (Figs. 1b and 2). The spatial and temporal thresholds were chosen so as they were the minimum values necessary to retain a large enough sample size for further analyses. Using the spatial and temporal thresholds described, only Clade 1 had enough spatially and temporally associated isolates to warrant further analyses.

Clustering of inter-isolate genetic distances

The available data for the isolates-REA type, sampling location (district where the sampling took place) and sampled host - were used to define groups of isolates and the within- and between-group genetic distances were examined to determine whether there was an association. The concatenated sequence of VPs of each isolate was compared to one another to generate an inter-isolate genetic distance (using the p-distance-defined as the proportion of the sites that differ between two sequences).

The observed difference between the mean intra- and inter-group genetic distances was calculated where the groups were defined, separately, by host species sampled, sampling location and REA type. To determine whether each observed difference could have arisen by chance alone, the isolate data was shuffled and the difference re-calculated. The shuffling was repeated 10,000 times to generate null distributions of observed differences. The associations were considered significant if the observed metric fell outside the lower (2.5%) and upper (97.5%) quantiles of the null distribution. Importantly, any species signature is likely to be nested within a spatial one, since regional localisation of bTB is known. To account for this, only comparisons that were between isolates sampled in the same district were included in the clustering analyses using the host species sampled.

Phylogenetic Analyses

The Bayesian Evolutionary Analysis by Sampling Trees (BEAST v1.8.4 - Drummond & Rambaut [36]) software was used for a phylogenetic analysis of the isolates’ sequences combined with their sampling years. BEAST was used to estimate the phylogenetic tree topology, substitution rate and date of the Most Recent Common Ancestor (MRCA) for the sampled M. bovis population. A BEAST analysis requires the existence of a clock-like substitution process. Additional analyses, as conducted by Firth et al. [37], were used to examine whether a clock-like process could have produced the inter-isolate variation (See Additional file 1: Supplementary text 1.3).

Models selected in a BEAST analysis may significantly impact the results. Care must be taken to select appropriate models for the substitution process [38] and the underlying population dynamics [39]. A series of BEAST analyses were completed in a hierarchical fashion to explore the different models available; for each analysis a chain length of 500,000,000 steps, sampled every 50,000 steps, was used and three replicates were completed. Following the removal of a 10% burn-in, the posterior distributions were examined to determine which structure of BEAST analysis best described the isolate data. Different analyses were compared based upon the log likelihood scores, model convergence and posterior support of parameters (assessed using TRACER v1.6 [40]), path sampling and stepping stone analyses (See Additional file 1: Supplementary text 1.4). In addition, the biological feasibility of the results was examined for each analysis.

The selected BEAST analysis used the Hasegawa-Kishino-Yano (HKY) substitution model, a relaxed clock model, drawing from an exponential distribution, and the Gaussian Markov Random Field (GMRF) Bayesian Skygrid population model. The HKY substitution model allows variable base frequencies, transition and transversion rates to be estimated [41]. A relaxed clock model enabled the estimated substitution rate to vary across the branches of the phylogenetic tree; the extent of this variation was modelled using an exponential distribution. The GMRF Skygrid model is a flexible model that is able to estimate changing population dynamics over the course of a phylogenetic history [39]. In a BEAST analysis, population dynamics are estimated based on the structure of the phylogenetic tree according to coalescent theory [42].

An additional Discrete Ancestral Trait Mapping (DATM) analysis [38, 43] for two states was implemented in the BEAST analysis. According to the host species sampled, isolates were assigned either a cattle or wildlife state. Based upon the states of the tips of the phylogenetic tree (the isolates) the DATM estimates the ancestral states in the phylogeny, and as such the most likely sources of infection within the sampled M. bovis population. A comparison was made between a symmetric and asymmetric DATM analyses in BEAST using the spatially and temporally matched isolates. The former symmetric analysis refers to the state transition matrix being symmetric; this analysis estimates a single parameter (in a two state analysis), the transmission rate of the pathogen from one state to another. The asymmetric analysis has two inter-state transmission parameters and as such can be used to determine whether there is a directional bias in the exchange; is the pathogen jumping from one population into another more often than in the other direction?

The influence of the selection of prior distributions for the parameters estimated in the BEAST analyses, described above, was investigated by running an analysis where the data were removed and only the prior distributions sampled. It was shown that the selected prior distributions were conservative and that the data provided a strong signal for the parameter estimations of our model (See Additional file 1: Supplementary text 1.5).

Structure in the sampled M. bovis population

There were four recognisably distinct clades formed by the 321 isolates sampled in New Zealand that were regionally localised (Fig. 1). A total of 3449 VPs were found. Long distance translocation and establishment of new foci of infection was evident, when the genetic structure of the population was considered. Clade 2 (Fig. 1 – red), although mostly found in New Zealand’s North Island, has an established foci of infection involving both cattle and wildlife on the South Island. Clade 1 (Fig. 1 – blue) isolates were mostly situated on the South Island of New Zealand, providing a densely sampled genetically similar set from which to select the spatially and temporally associated isolates for further analysis (Figs. 1 and 2). Clade 3 (Fig. 1 – gold) included eight wildlife and five cattle isolates, found across a broad spatial range in the southwest of New Zealand’s South Island. Clade 4 (Fig. 1 – green) also included thirteen cattle and three wildlife isolates, which were sampled from two locations <20 km apart in the south of New Zealand’s South Island.

Clustering of inter-isolate genetic distances

The inter-isolate genetic distance distribution of the spatially and temporally associated isolates from clade 1 was examined. Isolates of the same REA type were, on average, more genetically similar than those of different types. This difference was reflected in lower average within- than between-group genetic distances, when groups were defined by REA types (Fig. 3b). In addition, diversity was evident in the within-group distances demonstrating the added resolution of WGS data. The observed difference between the mean inter- and intra-group genetic distances was unlikely to have arisen by chance when the isolates were grouped by their REA type, sampling location or sampled species (Fig. 3b and c). In contrast to the lower within- than between-group genetic distances observed when groups were defined by REA type or sampling location (as was demonstrated by the positive observed difference (Fig. 3b and c)), when groups were defined by the host species sampled, the within-group distances were higher than the between-group distances (Fig. 3d). These higher within-group distances resulted in a negative observed difference, which was unlikely to have arisen by chance as it fell just outside the 95% bounds of the generated null distribution. This negative difference may be caused by lower within-outbreak distances resulting from sampling local outbreaks (involving cattle and wildlife) that are separated in space.

Substitution rate estimation

Using a bootstrapping procedure the posterior distributions resulting from BEAST analyses incorporating different population models were compared (Fig. 4b). For each pairwise posterior comparison, a distribution of differences was generated by calculating the difference between single point estimates, that were sampled proportionately from each of the two posterior distributions. If similar distributions are compared using this pairwise comparison, the calculated differences between point estimates drawn randomly from each distribution will be close to zero. The paired posterior distributions were not significantly different; the distribution of calculated differences resulting from each pairwise comparison overlapped with zero. The Skygrid population model, which had a high likelihood in the model selection procedures (See Additional file 1: Supplementary text 1.4) and agreed well with the other population models used (Fig. 4), estimated the substitution rate of the sampled M. bovis population to be 0.53 (2.5% Lower: 0.22, 97.5% Upper: 0.94) events per genome per year.

Using the Skygrid population model, the MRCA to the sampled M. bovis population was estimated to have been circulating in 1859 (2.5% Lower: 1525, 97.5% Upper: 1936). Binney et al. [44] recently established that a large number of cattle were imported into New Zealand in the 1860s, mostly originating from Australia and the United Kingdom. The structure of clade 1 and of the full maximum likelihood phylogeny (Fig. 1a) aren’t indicative of a single introduction event into New Zealand. The sampling window used in this study was narrow, relative to the phylogenetic history of the sampled population. Via simulation it was shown that estimates of the substitution rate were robust to the shortening of the sampling window; estimates increasingly lacked precision but retained accuracy (data not shown).

Ancestral traits analysis

Different selections of isolates from clade 1 (all isolates, temporally and spatially matched isolates, or 30 random cattle and wildlife pairs) were used in separate DATM BEAST analyses. According to the path sampling likelihood values, the symmetric model (equal rates from cattle to wildlife and vice versa) was favoured for the matched isolates (Fig. 5a). The asymmetric (different rates) was favoured when the DATM analysis was based on all the clade 1 isolates or the randomly paired isolates (Fig. 5a). The DATM analyses were able to provide an estimate for the overall state transition rate (Fig. 5b). For the analyses based on all the clade 1 isolates or randomly matched isolates the asymmetric model provided directional estimates of the state transition rates (Fig. 5c and d). Using all isolates or the random cattle and wildlife pairs, a dominant direction of transmission from wildlife to cattle was estimated (Fig. 5c and d). Using the spatially and temporally matched isolates the symmetric model out-performed the asymmetric model and directional estimates weren’t available. The difference in the support for the symmetric versus asymmetric model was a result of the isolates selected and therefore shows a strong influence of sampling. Without knowing which sample set is most representative it is difficult to have confidence in the directional state transition rates.