Full genome SNP-based phylogenetic analysis reveals the origin and global spread of Brucella melitensis

Brucellosis is an important zoonotic disease that affects both humans and animals. We sequenced the full genome and characterised the genetic diversity of two Brucella melitensis isolates from Malaysia and the Philippines. In addition, we performed a comparative whole-genome single nucleotide polymorphism (SNP) analysis of B. melitensis strains collected from around the world, to investigate the potential origin and the history of the global spread of B. melitensis. Single sequencing runs of each genome resulted in draft genome sequences of MY1483/09 and Phil1136/12, which covered 99.85% and 99.92% of the complete genome sequences, respectively. The B. melitensis genome sequences, and two B. abortus strains used as the outgroup strains, yielded a total of 13,728 SNP sites. Phylogenetic analysis using whole-genome SNPs and geographical distribution of the isolates revealed spatial clustering of the B. melitensis isolates into five genotypes, I, II, III, IV and V. The Mediterranean strains, identified as genotype I, occupied the basal node of the phylogenetic tree, suggesting that B. melitensis may have originated from the Mediterranean regions. All of the Asian B. melitensis strains clustered into genotype II with the SEA strains, including the two isolates sequenced in this study, forming a distinct clade denoted here as genotype IId. Genotypes III, IV and V of B. melitensis demonstrated a restricted geographical distribution, with genotype III representing the African lineage, genotype IV representing the European lineage and genotype V representing the American lineage. We showed that SNPs retrieved from the B. melitensis draft full genomes were sufficient to resolve the interspecies relationships between B. melitensis strains and to discriminate between the vaccine and endemic strains. Phylogeographic reconstruction of the history of B. melitensis global spread at a finer scale by using whole-genome SNP analyses supported the origin of all B. melitensis strains from the Mediterranean region. The possible global distribution of B. melitensis following the ancient trade routes was also consistent with whole-genome SNP phylogeny. The whole genome SNP phylogenetics analysis, hence is a powerful tool for intraspecies discrimination of closely related species.


Background
Brucella spp. are listed as a category B potential biological warfare agent [1]. A low infectious dose of 100 to 1,000 organisms is sufficient to cause an infection. The mechanism of transmission, through aerosols or food chains, makes them easily transmissible to both humans and animals. Extensive studies have been done in the past exploring the potential of Brucella spp. as candidate biological weapon agents [2]. A biological warfare attack model built in 1997, based on an aerosol attack using B. melitensis in an urban population, was estimated could cause an economic loss of $477.7 million per 100,000 persons exposed [3]. Owing to the potentially insidious nature of a bioterrorism attack, the intentional release of any biological warfare agent could be difficult to detect. Tracking unusual strains in a particular region could serve as the first line of defence. Rapid detection of potential bioterrorism activities that usually starts with few indistinguishable cases from a natural outbreak, however, is difficult. It relies heavily on precise and sensitive methods, as well as the availability of information about the organism worldwide. Sporadic outbreaks of brucellosis occur worldwide, the probable geographical origin of the source of infection, however, is difficult to determine.
Current approaches to study B. melitensis genome are often limited to specifically targeted genes, for example, multi-locus variable-number tandem repeat (MLVA) analysis, multi-locus sequencing (MLS) and single nucleotide polymorphisms (SNPs) analysis are insensitive to detect new mutations from new clades and therefore, are not very useful [4,5]. Complete genome sequencing is the ideal. The methods for achieving it, however, were previously laborious, expensive and usually required time-consuming gap closing of genome assemblies. It was therefore only used for in-depth evolutionary analysis or screening of potentially useful SNPs or loci to develop assays such as the MLVA, MLS, multiple real time and gene sequencing [6]. With the development of next-generation sequencing (NGS) technologies, full bacterial genome sequencing has become highly accessible [4,[7][8][9][10][11]. High quality draft bacterial genomes can be rapidly obtained. Unlike the taxonomically informative or canonical SNP-based approaches, whole-genome sequencing serves as a robust and unbiased method to resolve intraspecies relationships for closely related species such as Brucella spp. [12,13] and Bacillus anthracis [14,15].
To date, there are only a limited number of full genome sequences of B. melitensis isolates from different geographical locations available in publicly accessible databases. Only one draft full genome sequence of B. melitensis isolated from the Southeast Asia (SEA) region is available [16]. We identified and sequenced two human B. melitensis strains, MY1483/09 and Phil1136/12, which originated from Malaysia and the Philippines, respectively. We demonstrate here that a SNP-based phylogenetic analysis of the draft genome sequence obtained from a single sequencing run was sufficient to resolve the intraspecies genetic relationships of geographically distinct B. melitensis strains. In addition, whole-genome-based molecular epidemiology tracking allows the elucidation of the potential origins and spread of B. melitensis.

General features
The full genome sequences of the Brucella isolates were generated from a single sequencing run using the 316 chips of the Life Technologies Ion Torrent PGM platform. The resulting reads were de novo assembled into two unclosed draft genomes, comprising 28  A total of 3,355 and 3,408 protein encoding genes were predicted and allocated in the chromosomes of the MY1483/09 and Phil1136/12 isolates, respectively, using GeneMarkS v4.10. The predicted functions of 2,920 genes of MY1483/09 (87% of total predicted genes) and 2,978 genes of Phil1136/12 (87% of total predicted genes) were obtained by searching against the NR, Swiss-Prot and KEGG databases. Three copies of rRNA in both isolates were identified using rRNAmmer v1.2. A total of 50 and 51 copies of tRNA were predicted using tRNAscan v1.23 in MY1483/09 and Phil1136/12, respectively ( Figure 1A and B).

Phylogenetic analysis
The whole-genome-based SNP phylogenetic analysis was used to infer the relationships between B. melitensis strains collected from diverse geographical locations. We compared the draft genome sequence of the two isolates obtained in this study against 49 strains of B. melitensis (complete genome and draft genome) available in online databases at the time of the analysis. B. melitensis 16 M was used as the reference genome for sequence dataset, containing only the assembled sequences for the SNP discovery using SNPsFinder [17]. As a previous study indicated a close phylogenetic relationship between B. melitensis and B. abortus [4], the genomes of two B. abortus strains were selected as the outgroup in this analysis. The paralogous genes were removed from the analysis. A comparison of all Brucella strains, including the outgroup strains, with a mismatch cut-off of 8 bp, yielded a total of 13,728 polymorphic nucleotide sites (Additional file 1: Table S1). A phylogenetic tree was reconstructed using the identified SNPs shared among these genomes. The phylogenetic tree, together with the details of the isolates, including the geographical information, the host and the biovar were shown in Figure 2A. All of the B. melitensis strains, except 16M13W and S66, segregated into five major genotypes that corresponded spatially with the potential origins of the isolates. B. melitensis Ether (biovar 3) was basal to the B. melitensis clade. B. melitensis 16 M (biovar 1) and 63/9 (biovar 2) were separated after divergence from the Ether strain, demonstrating a similar tree topology that previously reported [4]. Incorporating the tree topology and the age of the B. melitensis lineage previously reported [4], we denoted B. melitensis Ether and its related strains as genotype I, B. melitensis 63/9 and its related strains as genotype II and B. melitensis 16 M and its related strain as genotype V. The B. melitensis genotypes III and IV were new genotypes that had not been reported, consisting mainly of recent isolates. Their genetic relationships with other B. melitensis strains have therefore not been well described [16]. The overall geographical distribution of the B. melitensis genotypes and clades is illustrated in Figure 2B.
Genotype I represented the most basal lineage of the B. melitensis strains. It comprised four out of the 49 B. melitesis strains used in the study. Based on the geographical origin of the strains, most of the genotype I B. melitensis strains were isolated from Italy (Ether, F15/06-7 and F5/07-239A) and one isolate, UK31/99 was collected from Egypt. The genotype I strains contained a number of SNPs, ranging from 2,618 to 2,762 in comparison to the B. melitensis 16 M. By conducting the SNP analysis, we identified 936 polymorphic sites that distinguished genotype I strains from other genotypes. A pairwise comparison of the genotype I isolates revealed the presence of 745 intragenotypic SNPs.
Genotype II comprised 30 out of the 49 B. melitensis strains used in the study. It represented the largest B. melitensis genotype, with isolates collected from diverse geographical locations, including African, European and Asian countries. It exhibited a ladder-like phylogram, suggesting a possible single introduction of genotype II strains into the region. The genotype II strains possessed SNP differences ranging from 2,187 to 2,521 in comparison to B. melitensis 16 M. Within genotype II, the isolates further segregated into several clades. A small clade, denoted group a, consisted of four B. melitensis strains that emerged from B. melitensis UK3/06 (Cyprus). The genotype IIb strains were collected from distantly related geographical locations, with two of the strains having been collected from European regions, Kosovo and the United Kingdom (UK), one strain from Turkey (Asia) and one strain from Somalia (Africa). A pairwise comparison of the genomes revealed the presence of 295 SNPs within the IIb strains. We identified 67 SNPs that were specific to B. melitensis genotype IIb.
UK29/05, a B. melitensis strain collected from the UK, emerged from genotype IIb. It occupied the internal node on the phylogenetic tree prior to the emergence of   Genotypes III, IV and V emerged from the same common ancestral strains and demonstrated strong geographical clustering. Genotype III consisted of strains belonging to the African group, with the exception of UK22/04, which was collected from the UK. The isolates from genotype III had SNPs that numbered from 1,637 to 1,698 in comparison to B. melitensis 16 M. Based on the whole-genome SNP analysis, we identified 255 SNPs that were specific to genotype III strains. A pairwise comparison of genotype III strains revealed the presence of 1,388 intragenotypic SNPs. Genotype IV comprised isolates collected from Norway, Portugal and Malta (European clade). We identified 97 clade-specific SNPs that were shared among the genotype IV strains. A pairwise comparison of the genotype IV isolates showed 163 intragenotypic SNPs. Genotype V comprised the American strains collected from Argentina (CNGB1075, CNGB1120 and CNGB290), Mexico (Rev-1) and the USA (16 M strains except China 16 M). Within genotype V, the strains segregated into two groups, representing those from South America (Argentina) and those from North America (Mexico and USA). We identified 28 specific SNPs for the North American clades. However, we did not observe any specific SNPs for the South American clades.
Our analysis also included genome sequences of all of the vaccine strains of B. melitensis used worldwide. Rev-1 is a worldwide vaccine strain [18,19], whereas M5 and M111 are the vaccine strains currently used in China [8,10]. Rev-1 belonged to genotype V, grouping closely with 16 M and other isolates collected from Argentina. The M5 and M111 vaccine strains clustered with genotype II group f, the China clades. The M5 strain grouped closely with its parental strain, M28.

Discussion
Anthropogenic activities allow infectious agents to adapt and be disseminated more efficiently to diverse geographical locations by their host. By tracking specific genome modifications, inference of the evolutionary history of the B. melitensis strains can be made. The Brucella spp. genome comprises two chromosomes without any additional extrachromosomal genetic elements. The horizontal transfer of genes, a common mechanism contributing to bacterial genome diversity, is limited in Brucella spp., perhaps due to its strict intracellular niche in the host [11]. Genome modifications, including mutations and genes rearrangements, in Brucella spp. therefore, are the exclusive heritable elements during the evolution of the bacteria [20]. The accumulation of mutations or acquisition of genome specific structural changes reflects the microevolution events that drive the evolutionary process, through genetic drift or bacterial adaption, upon introduction into a new geographical region or host. Among the Brucella spp., B. melitensis is the most virulent species for humans. Here, by incorporating a wholegenome SNP-based phylogenetic analysis and identifying the geographical origins of the isolates, we demonstrated a strong phylogenetic link in B. melitensis in relation to its distinctive geographical origins.
The presence of Brucella spp. in the ancient Herculaneum [21] and the high prevalence (17.4%) of brucellosis in Herculaneum adults [22] demonstrates the prolonged history and continuous circulation of Brucella spp. in the Mediterranean area. A high diversity of B. melitensis strains were collected from countries at the border of the Mediterranean Sea [23][24][25] and their basal location (genotype 1) on the phylogram suggests that B. melitensis could have originated from this region. This is further supported by the recent discovery of a 700 year-old B. melitensis genotype I strain from Sardinia [26]. Ancient and Old World transcontinental trading was an early form of globalisation that allowed the exchange of goods among the continents. The infectious agents that attached to humans or goods, such as livestock and its derivatives spread in a parallel manner [27]. Based on the phylogenetic tree drawn in this study using whole-genome SNPs, it was shown that the B. melitensis genotype I gave rise to two major branches, representing B. melitensis strains collected from the East (genotype II) and West (genotype III/IV/V). The strains were likely to have been distributed along the major ancient trade routes, e.g., the ancient Silk Road, the Trans-Eurasia exchange and the Columbian exchange. The transcontinental spread of B. melitensis could therefore have marked the beginning of the diversification of the B. melitensis strains. The global spread of B. melitensis may have occurred following these ancient trading routes, perhaps through infected goat and sheep or their milk derivatives.
B. melitensis genotype II had the most diverse geographical distribution after diverging from genotype I. Following the phylogenetic ladder pattern, the distribution of genotype II strains expanded from West (Cyprus) to East (China). Its ladder-like tree topology suggests that the genotype II strains descended from a single common ancestor that originated from the Mediterranean area (Cyprus origin?) to form the Asian clades, as almost all Asian strains of B. melitensis belonged to this group. It is noted that B. melitensis BG2 (S27), isolated from Pakistan, occupied the ancestral node of group f (China clade), suggesting that the B. melitensis strain currently circulating in China may has a Pakistani origin.
Genotypes III, IV and V of B. melitensis demonstrated a more restricted distribution. They shared a common ancestor, which most likely diversified and spread towards the southern region of the Mediterranean, and formed genotype III (African lineage). The strains that spread northwards could have further diversified and adapted geographically forming genotype IV (European lineage) and genotype V (American lineage). Compared with the American region, genotype IV of B. melitensis strains were more geographically diverse (Portugal, Malta and Norway). Unlike Portugal and Malta, Norway is officially a brucellosis-free country [28], with no B. melitensis having been detected in sheep or goat herds to date [29]. Therefore, the Norway strain, F3/02, most likely represents an imported case that probably originated in the northern part of the Mediterranean region. There were limited numbers of sequences available for analysis from this region. Thus, we cannot rule out the sampling bias that may have affected the assignment of genotype IV to the northern Mediterranean area.
The genotype V strains were most probably introduced from Europe into the American continent by infected animals several hundred years ago, where B. melitensis could have continued to diversify into two geographically specific sub-lineages, genotype Va (North American lineage) and Vb (South American lineage). Nevertheless, only two strains (three genome sequences) and three strains were used to designate genotypes Va and Vb, respectively. The two genotype Va strains, 16 M (USA origin) [30] and Rev-1 (Mexico origin; vaccine strain developed in USA laboratories) [31,32], however, may not reflect the currently circulating B. melitensis strain in the North American region. Similarly, all three B. melitensis strains used for the designation of genotype Vb were isolated from Argentina. Although the origin of the Argentinian strains in this study (native or imported) remains inconclusive, the continuing circulation of B. melitensis in Mexico and Argentina, where B. melitensis is the principle cause of brucellosis [33], suggests the existence of native B. melitensis in the American continent. With the current data, we demonstrated that genotype V of B. melitensis was only detected in the American region, supporting the geographical specificity of genotype V. More field B. melitensis strains from diverse geographical regions of the American continent, however, are needed to conclude the designation of genotype V.
The SEA region is believed to have been free of B. melitensis until recently, similar to Australia, New Zealand, Northern Europe, the US and Canada [34,35]. Recent epidemiological studies of livestock and B. melitensis associated with human brucellosis, however, have suggested that the bacteria could have been transmitted and maintained unknowingly in the goat and sheep populations [36][37][38][39][40][41]. These SEA B. melitensis strains have not been well studied. Genetic markers, such as IS711 used for species characterisation in previous studies, are not sufficient to delineate intraspecies relationships of B. melitensis [40,41]. Here, using a comprehensive whole-genome-based SNP approach, the two B. melitensis strains obtained in this study clustered closely with the Thailand B. melitensis strains (F10/06-16) and formed a new SEA clade. Deciphering the origin and source of the SEA B. melitensis strains is highly challenging, owing to the existence of both local and imported goat breeds in SEA countries. We identified 110 SEA clade-specific SNPs that differentiated the SEA clade from other global B. melitensis strains, suggesting a common ancestor of the SEA B. melitensis strains in this region. A previous study suggested that the persistence of B. melitensis in wildlife could serve as a reservoir for endemic outbreaks [42]. The high genetic similarity (15 polymorphic SNPs) of the three SEA B. melitensis strains, however, suggests that the spread of SEA strains was relatively recent. It is therefore unlikely that brucellosis was caused by autochthonous transmission from local breeds in the wild. Moreover, Malaysia, the Philippines and Thailand are not major exporters of goats and sheep in the region. Therefore, it is plausible that the B. melitensis in this region was introduced from a common source through the importation of goats and sheep from elsewhere. The bacterial strains became established and adapted to the local environment, leading to the spread of B. melitensis to the local goat and sheep populations.
In view of the persistent infection caused by the vaccine strains in animals [43,44], their possible involvement in human infection is worth investigating. Data concerning the endemic strains and the use of vaccines in local animal husbandry, however, have not been well documented. It is possible that the vaccine strains used in animals could intentionally be used to cause harm to humans. There are several veterinary vaccine Brucella spp. strains currently in use. We showed here, by using a wholegenome SNP-based phylogenetic tree, that we were able to resolve and discriminate between the vaccine strains and the endemic strains. This is particularly important to help identify possible sources of infection and to monitor the potential misuse of live attenuated vaccines.
Through the whole-genome SNP analysis, two previously assigned B. melitensis strains, S66 and 16M13W, segregated into a group that did not belong to B. melitensis. S66 was previously reported as a B. melitensis sequence type 8 strain [45], while 16M13W is the 13-week mice-adapted strain of China 16 M B. melitensis [9]. Clustering of the 16 M13W strain into a different group from its parental strain suggests that it may not have originated from China 16 M. Although both China 16 M and 16M1W were identified as B. melitensis strains, they were not grouped with other previously designated 16 M strains (bv1 str. 16 M and bv1 str. 16 M_1). In contrast, they grouped with other B. melitensis strains collected from China (genotype II), raising concern that the China 16 M may not be the same 16 M strain used worldwide.
To date, serological tests, such as the agglutination test, and isolation of bacteria from blood cultures remain the gold standard for the diagnosis of human brucellosis. The interpretation of the results from such assays, however can be difficult due to the serological cross-reaction of Brucella spp. with other bacteria [46,47]. In addition, differential identification of closely related Brucella spp. is not possible using a single serological test. The PCR-based genome sequence amplification method is a good alternative, as it allows for rapid identification of the bacteria. Differential identification of highly similar Brucella spp., however is only achievable with multiple gene sequence analyses, such as MLS and MLVA. With the limited candidate genes used for MLS and MLVA typing, the intraspecies relationships of closely related species of B. melitensis are difficult to resolve. In comparison to these subtyping approaches, the whole-genome SNP-based phylogenetic analysis presented here provides better resolution power to resolve the genetic relationships between the B. melitensis species.

Conclusion
In summary, we have demonstrated that SNPs retrieved from the B. melitensis draft genomes were sufficient to resolve the interspecies relationships within B. melitensis strains. Our findings obtained using the genome SNP analyses revealed potential B. melitensis endemicity to the SEA region. The B. melitensis phylogeny presented here describes the relationships between the global B. melitensis strains and their evolutionary history. The whole-genome SNP-based approach provides an effective means of determining the native geographical origin of the pathogen and could serve as the basis for the reconstruction of the history of the global spread of B. melitensis at a finer scale. In addition, we have discriminated between the vaccine and endemic strains using the current whole-genome SNP approach. This will allow the rapid identification of pathogens and sensitive epidemiological follow-up in cases of release with intention to cause harm.

Genomic DNA sample preparation
The bacteria were cultured in Brucella broth and incubated at 37°C for 72 hours. Genomic DNA was extracted from 5 mL of culture using the QIAamp DNA Mini Kit (Qiagen, Germany). All sample handling was performed in a BSL-3 biocontainment laboratory at the Tropical Infectious Diseases Research and Educational Centre (TIDREC), University of Malaya. The genomic DNA concentration was quantified using the Qubit® dsDNA BR Assay Kit (Invitrogen, Life Technologies, USA).

Full genome sequencing
The genome library preparation was performed by using the Ion Xpress TM Plus Fragment Library Kit (Life Technologies, USA) following the manufacturer's protocol, 4471989 Rev E. The sheared DNA was separated using a 2% agarose gel (Promega, USA) and fragments that corresponded to 200 bp of a sequencing run (approximately 330 bp) were manually excised. The DNA was purified using the QIAquick Gel Extraction Kit (Qiagen, Germany). The library quality and concentration were assessed and determined using the 2100 Bioanalyzer (Agilent, USA) and High Sensitivity DNA Kit (Agilent, USA). The sequencing template was prepared using the Ion OneTouch TM 200 Template Kit V2 DL (Life Technologies, USA) according to the manufacturer's protocol. The Ion Sphere Particles were enriched using Dynabeads MyOne Streptavidin C1 beads (Invitrogen, Life Technologies, USA). The efficiency of the enrichment process was assessed using the Ion Sphere Quality Control Kit (Life Technologies, USA). Genome sequencing was undertaken using the Ion Torrent PGM sequencer (Life Technologies, USA). Outgroup Brucella abortus S19 Brucella abortus S19 GCA_000018725.1