Species identification and molecular typing of human Brucella isolates from Kuwait

Brucellosis is a zoonotic disease of major concern in Kuwait and the Middle East. Human brucellosis can be caused by several Brucella species with varying degree of pathogenesis, and relapses are common after apparently successful therapy. The classical biochemical methods for identification of Brucella are time-consuming, cumbersome, and provide information limited to the species level only. In contrast, molecular methods are rapid and provide differentiation at intra-species level. In this study, four molecular methods [16S rRNA gene sequencing, real-time PCR, enterobacterial repetitive intergenic consensus (ERIC)-PCR and multilocus variable-number tandem-repeat analysis (MLVA)-8, MLVA-11 and MLVA-16 were evaluated for the identification and typing of 75 strains of Brucella isolated in Kuwait. 16S rRNA gene sequencing of all isolates showed 90–99% sequence identity with B. melitensis and real-time PCR with genus- and species- specific primers identified all isolates as B. melitensis. The results of ERIC-PCR suggested the existence of 75 ERIC genotypes of B. melitensis with a discriminatory index of 0.997. Cluster classification of these genotypes divided them into two clusters, A and B, diverging at ~25%. The maximum number of genotypes (n = 51) were found in cluster B5. MLVA-8 analysis identified all isolates as B. melitensis, and MLVA-8, MLVA-11 and MLVA-16 typing divided the isolates into 10, 32 and 71 MLVA types, respectively. Furthermore, the combined minimum spanning tree analysis demonstrated that, compared to MLVA types discovered all over the world, the Kuwaiti isolates were a distinct group of MLVA-11 and MLVA-16 types in the East Mediterranean Region.


Introduction
Human brucellosis, a common zoonotic disease, is a neglected, under-recognized infection of widespread geographic distribution and globally about 500,000 cases occur annually [1]. The highest incidence of human brucellosis is recorded in the Middle East and Central Asia [2,3]. It is among the most commonly reported infectious diseases in Kuwait  In this study, we have identified the predominant Brucella spp. infecting humans in Kuwait by 16S rRNA gene sequencing and real-time PCR HRM analysis, and genotyped the isolates using ERIC-PCR and MLVA techniques. The genetic relatedness of Kuwaiti strains with the strains isolated internationally was determined by comparing the MLVA types using minimum spanning tree (MST) analysis.

Materials and methods
Brucella strains and DNA isolation A total of 75 Brucella strains (BRU001-BRU118) isolated from 75 patients and cultured on plates at the Clinical Microbiology Laboratories of Infectious Diseases, Mubarak Al-Kabeer, Farwaniya and Amiri Hospitals in Kuwait. A loopful of bacterial colonies from each plate was suspended into 1 ml phosphate buffered saline (pH 7.0) and heated at 95˚C for 10 min in a water bath. The genomic DNA was isolated from the heated specimens using the QIAamp DNA Mini Kit (Qiagen, Valencia, CA) according to the manufacturer's instructions. The quantities and purities of isolated DNA were determined using an Epoch Spectrophotometer (Biotek, Winooski, VT) and Qubit Fluorometer (Qubit dsDNA BR Assay Kit, Life Technologies, Carlsbad, CA). The isolated DNA was stored at −80˚C until further use.

Amplification and sequencing the 16S rRNA gene
A 500 bp region of 16S rRNA gene from Brucella genomic DNA was amplified using the MicroSeq 1 500 16S rDNA PCR kit (Applied Biosystems, Grand Island, NY) according to the manufacturer's instructions. In brief, the PCR reaction mixtures contained 15 μl of genomic DNA (25 ng) and 15 μl of 2x master mix (the universal primers, AmpliTaq 1 gold DNA Polymerase, Buffer, MgCl2 and dNTP mix) from the kit, and the PCR was performed in a Gen-eAmp 1 PCR System 9700 (Applied Biosystems, Grand Island, NY). PCR mixtures containing positive (DNA isolated from Escherichia coli) and negative controls (nuclease-free water), included in the kit, were also prepared. The conditions used for target amplification were as follows: 95˚C for 10 min, 30 cycles of 95˚C for 30 s, 60˚C for 30 s, and 72˚C for 45 s, and a final extension step at 72˚C for 10 min. The PCR products were visualized on a DNA 1000 gel in a Bioanalyzer (Agilent 2100, Santa Clara, CA). The amplified products were purified by adding 2μl of ExoSAP-IT 1 to 5μl of the PCR product and incubated for 15 min at 37˚C, followed by heat inactivation of the enzyme at 80˚C for 15 min. Cycle sequencing was performed with the purified PCR products using the MicroSeq 1 500 sequencing kit (BigDye 1 Terminator v1.1 chemistry) as per the kit protocol. The thermal cycler was programmed at 96˚C for 10 s, 50˚C for 5 s and 60˚C for 4 min for 25 cycles. The extension products were purified with Performa 1 DTR Gel Filtration Cartridges (Edge Biosystems, Gaithersburg, MD) and sequenced on an ABI 3130xl automated Genetic Analyzer (Applied Biosystems, Foster City, CA). The ABI files were opened with the Sequencing Analysis software (Applied Biosystems, Foster City, CA) for the quality assessment. DNA sequences in fasta format were further uploaded in the BioNumerics version 7.5 software (Applied Maths, Sint-Martens-Latem, Belgium) and submitted to the Ribosomal Database Project (RDP) [42] via the RDP plugin for similarity scores and percent identity calculation. The DNA sequence data have been submitted to the GenBank database under accession numbers MF164063 to MF164137.

Real-time PCR and HRM analysis
The real-time PCR assays were performed using Brucella genus-specific forward (O1: 5'-TC CGCAAGCTTCAAGCCTTCTATC-3') and reverse (O2: 5'-GGCGTGTCTGCATTCAACGTA ACC-3') primers [43], and B. melitensis-specific forward (BF: 5'-CATGCGCTATGTCTGG TTAC-3') and reverse (BMR: 5'-AGTGTTTCGGCTCAGAATAAT-3') primers [44]. The real-time PCR mixture was prepared by adding 5 μl of 10× PCR Buffer II, 0.7 μM final concentration of the forward and reverse primers (1μl each), 10 ng of DNA template (2 μl), SYBR Green I (2 μl), MgCl2 (2.4 μl) and nuclease free water (Qiagen, Germany) to a total volume of 25 μl per reaction. The real-time PCR was performed in a Light Cycler 1 2.0 (Roche Diagnostics GmbH, Mannheim, Germany) with an initial denaturation step of 95˚C for 10 min, followed by 35 cycles of 95˚C for 15 s and 65˚C for 10 s, with acquisition of data at 72˚C for 15 s in the green channel (excitation at 470 nm and detection at 510 nm). After amplification, the HRM analysis was performed between 65˚C and 95˚C at the rate of 0.1˚C.

ERIC-PCR
ERIC-PCR was performed with the primers ERIC1R (5 0 -ATGTAAGCTCCTGGGGATTCAC-3') and ERIC2 (5 0 -AAGTAAGTGACTGGGGTGAGCG-3') as described previously [36]. In brief, the PCR reaction mixtures contained 6 μl of 5x HOT FIREPol 1 Blend Master mix (Solis BioDyne, Estonia), 10 ng DNA, 0.15 μM of each primer in a total volume of 30 μl. PCR conditions were 95˚C for 12 min; 45 cycles of 95˚C for 45 sec, 35˚C for 1 min and 70˚C for 10 min; and a final step at 70˚C for 20 min. The amplification products were visualized on a DNA 7500 gel by the Bioanlayzer (Agilent 2100) and the band patterns were analyzed by the BioNumerics version 7.5 software (Applied Maths). In-house validation of the assay was done for the reproducibility employing three technical and biological replicates. Clustering was performed using the Dice similarity coefficient (optimization of 0.5%, tolerance of 1% and active zones of 10-78%) and unweighted-pair group method using arithmetic averages (UPGMA). Clusters were further classified based on ca. 80% of similarity [36]. The Discriminatory Index (DI) was calculated by the online tool (insilico.ehu.es/mini_tools/discriminatory_power/index.php).

MLVA typing
Amplification of the 16 VNTR loci for MLVA typing was performed, using primer pairs specific for each locus, according to the methods described previously [27][28][29][30][31][32][33][34]. In brief, PCR was performed in a total volume of 20 μl containing 1 ng of DNA, 1× PCR Master Mix (Solis Bio-Dyne, Tartu, Estonia), and 0.6 μM of each forward and reverse primers. Thermal cycling was performed in a GeneAmp 9700 Thermal Cycler (Applied Biosystems) by initial heating at 95˚C for 3 min, followed by 30 cycles of denaturation at 95˚C for 30 s, annealing at 60˚C for 30 s and extension at 72˚C for 50 s. The PCR products were visualized on a DNA 1000 gel by the Bioanlayzer (Agilent 2100) for accurate band size estimation. DNA 1000 ladder (100-1000 bp fragment size) from the Agilent kit was used as a control to estimate the band sizes of PCR products. Band sizes were converted into number of tandem repeat units according to the 2013 Brucella allele assignment table (version 3.6, available at http://mlva.u-psud.fr). The data set was submitted to the Brucella MLVA Database (http://mlva.u-psud.fr) for genotype identification on the basis of Panel 1 (MLVA8-type), Panel 1+2A (MLVA11-type) and Panel 1+2A +2B (MLVA16-type). All isolates were identified at Brucella spp. level by the MLVA8 typing scheme. The repeat numbers were imported to BioNumerics 7.5 software as character sets for cluster analysis based on categorical coefficients and UPGMA on similar weight basis. The discriminatory power of MLVA markers were calculated by the Hunter and Gaston Diversity Index (HGDI) via the online tool V-DICE available at the HPA website (http://www. hpabioinformatics.org.uk/cgi-bin/DICI/DICI.pl).
In order to determine the genetic relatedness of B. melitensis strains of our study with the strains available in the MLVA database, the VNTR copy numbers of 827 MLVA types from four major clades of the world (Africa, Americas, and East and West Mediterranean Regions) were downloaded from the MLVA database. Neighbor-joining minimum spanning trees were constructed from the datasets of the downloaded strains and the strains from Kuwait using BioNumerics 7.5 software.

Results
16S rRNA gene sequencing and real-time PCR PCR of the 16S rRNA genes of all the 75 isolates showed bands at about 500 bp. The positive control containing E. coli DNA also produced a band at~500 bp, whereas the negative control did not produce any bands (S1 Fig). The nucleotide lengths of the PCR products obtained after 25 cycles of sequencing PCR were~500 bp long (S1 Table). The quality value of each base call was >20 and an overall specimen score ranging from 26-46 was obtained for all the samples through the Sequencing Analysis software. Submission of all the 75 sequences to RDP database via BioNumerics identified the isolates as B. melitensis, exhibiting similarity values of 0.827 to 1.000 and 90-99% sequence identity (S1 Table). The results of real-time PCR assays with genus-and specific-specific primers confirmed that that all isolates belonged to genus Brucella, and species B. melitensis (data not shown).

ERIC-PCR
The ERIC-PCR products were well resolved in the Agilent Bioanlayzer, yielding peaks at corresponding bands (S2 Fig). The ERIC primers generated polymorphic band patterns in all the 75 isolates with varying numbers (13 to 39) and sizes of bands (73 bp to 5000 bp). The in-house validation results indicated the reproducibility of the assay by producing bands of similar intensities and length in all the technical and biological replicates (S3 Fig). The band profiles produced a dendrogram of 75 branches (Fig 1), which suggested that each sample was of a unique type. The DI calculated based on this finding was 0.997 confirming the high discriminatory power of the technique. The cluster classification based on~80% similarity divided all the ERIC genotypes into two major clusters, A and B. Cluster A consisted of 9 ERIC genotypes (A1-A9) corresponding to 9 individual strains. Cluster B comprised of 13 sub clusters (B1-B13) amongst which B5 formed the largest cluster with 51 strains (68% of total population) sharing more than 80% similar ERIC profiles (Fig 1), followed by B3 consisting of 5 strains. Clusters B1 and B4 were composed of 3 strains and the remaining B2, B6, B7, B8 and B9 had only one strain each.

MLVA typing
PCR amplification products were obtained with the sixteen MLVA primers in the presence of DNA from all 75 Brucella strains. However, varying number of alleles were detected, which ranged from 11 alleles for Bruce04 to only one allele for Bruce45 (Table 1). Overall, MLVA primers exhibited high discriminatory power, as observed from the HGDI value of 0.942. However, the individual set of primers varied in HGDI values ranging from as low as 0.000 (Bruce45) to as high as 0.888 (Bruce04) ( Table 1).
The VNTR numbers derived through the banding pattern obtained from Panel 1 (MLVA-8) markers matched with B. melitensis in the MLVA database. Furthermore, the Panel 1 markers divided the 75 isolates into ten MLVA-8 genotypes, three of which have been reported previously and seven genotypes that were unique (Fig 2). Among the previously reported MLVA-8 genotypes, the genotype number 45 was the predominant type in Kuwait (40 strains), followed by 64 (22 strains) and 57 (one strain). Seven unique genotypes, not reported in the MLVA-8 database, were named as K1, K2, K4, K5, K6, K7 and K8 (Fig 2). Two strains each belonged to K1 and K5 genotypes, and three strains each in K2 and K8; whereas the remaining three unique genotypes had only one strain, i.e. K7-BRU084, K2-BRU115, K4-BRU006 (Fig 2).
In order to determine the genetic relatedness of B. melitensis strains of our study with the strains available in the MLVA database, the VNTR copy numbers of 827 MLVA types from four major clades of the world (Africa, Americas, and East and West Mediterranean Regions) were downloaded from the MLVA database (S2 Table). Neighbor joining minimum spanning trees were constructed from the datasets of the downloaded strains and the strains from Kuwait. This method suggested that the MLVA-8 types found in Kuwait were present as two major groups i.e. 45 and 64, and the former originated from the latter (Fig 5). Both the groups appeared as a branch in the East Mediterranean clade tree (Fig 5). Interestingly, the genotype 64 showed a minimum distance of 1.0 with MLVA-8 type-43 of a strain isolated in the United    Arab Emirates, a country in the East Mediterranean Region [41]. Amongst the new MLVA-8 types, five strains (belonging to K1, K2, K4 and K8) clustered with the East Mediterranean region, two strains (belonging to K5) clustered with the African clade and one strain (belonging to K6) formed a branch in the West Mediterranean Region (Fig 5). In the MST made on the basis of MLVA 11 (panel 1+2A) dataset, a clear cut branching of Type 45 and 64 into several genotypes was observed (Fig 6). The MLVA-16 genotypes appeared as a profusely divided branch of the East Mediterranean clade (Fig 7).

Discussion
In this study, we have reported the results of testing four molecular methods to identify and genotype Brucella spp. infecting humans in Kuwait. The sequencing of 16S rRNA gene suggested that all 75 isolates were B. melitensis. Sequencing of the 16S rRNA gene is widely used as a speedy and accurate tool for bacterial identification, including Brucella [17,[45][46][47]. The total length of the 16S rRNA gene is~1500 bp, but the bacterial genera and species have been identified based on the first 500 bp region [48]. Suitability of partial 16S rRNA gene sequencing has been demonstrated for identification of dangerous pathogens, including B. melitensis [48,49,50]. The MicroSeq 1 500 microbial identification system is a robust and accurate tool for this purpose. It saves time and resources used for full-length gene sequencing. A straight forward workflow allows easy and fast handling and very quick results on the same day. The BioNumerics software aided in maintaining a comprehensive database that could be directly linked to the RDP database and easily accessible Sab and Similarity scores could be obtained rapidly. Although, the inability for intraspecific differentiation of Brucella is a limitation of this  Table). The MST was constructed using BioNumerics 7.5 software. The strains were categorized on the basis of their geographic location and differentiated through color codes.  [19]. Similar work was carried out recently by Zahidi et al. in Malaysia [20]. They concluded that the methodology of real-time PCR combined with HRM analysis was a fast, accurate and cost effective for identification of B. melitensis [20]. Monitoring the results in real-time saves from additional steps of gel electrophoresis. Moreover, detection of B. melitensis as a prevalent species in Kuwait would further aid to design eradication strategies, like vaccine development.
The genetic discrimination of Brucella remains a challenging task owing to its genetic homogeneity. However, typing of B. melitensis isolates is highly desirable for contact tracing  Table). The MST was constructed using BioNumerics 7.5 software. The same color codes, as given in Fig 5, were used to differentiate between the strains isolated from Kuwait and other parts of the world.
https://doi.org/10.1371/journal.pone.0182111.g006 and epidemiological outbreak investigations in Kuwait and Middle East. ERIC-PCR has been used as a typing tool in the past for Brucella [35] and other bacterial species [36]. Our results suggest that ERIC-PCR was a reliable test for identifying genetic differences within B. melitensis spp. None of our strains shared cent per cent resemblance with each other, hence 75 ERIC genotypes were identified in the region. Owing to the high resolution power of the technique strains sharing only 100% similarity were counted as similar types [36]. The resolution of bands was further improved by the use of Agilent Bioanalyzer. In the present investigation band profiles with high quality sizing resolution (1000-7500 bp: 15%), sizing accuracy (± 10%  Table). The MST was constructed using BioNumerics 7.5 software. The same color codes, as given in Fig 5, were used to differentiate between the strains isolated from Kuwait and other parts of the world. https://doi.org/10.1371/journal.pone.0182111.g007 Molecular characterization of Brucella in Kuwait CV) and sizing reproducibility (5% CV) were obtained, which are essential requirements in fingerprinting based differentiation [51,52]. In our study, the Agilent Bioanalyzer's band profiles along with the most advanced version of BioNumerics 7.5 formed a model experimentanalysis combination. The application of DICE algorithm in the BioNumerics software scored the bands as present or absent and created dendrograms that could be used for cluster analysis. The dendrogram of B. melitensis was classified into two clusters, A and B, based on~80% similarity [36]. The maximum number of strains (n = 51; 68% of total population) were present in the cluster B5. The closeness among 68% of strains of Brucella could be attributed to the generation of a large number of identical bands corresponding to the conserved region. However, the typeability was calculated on the basis of limited number of polymorphic bands. Our findings confirm the previous reports of Mercier and co-workers [35] that ERIC-PCR is capable of differentiating between the Brucella strains on account of even limited number of polymorphic fragments. In a relatively recent study on the highly homogenous Corynebacterium pseudotuberculosis, the ERIC-PCR successfully typed the various strains with high power of discrimination and reproducibility [36].
The species identification and genotyping of all 75 Brucella isolates was further extended using the MLVA technique [53][54][55] and the online Brucella database (http://mlva.u-psud.fr/ brucella/). This database has been extensively used by other investigators for identification and typing of Brucella [27-33, 40-42, 56], and it is regularly updated. The latest version (released on May 16, 2016) contains data for more than 4000 Brucella strains of various species. In our study, all isolates were identified as B. melitensis by the MLVA-8 (Panel 1) typing scheme. However, for the purpose of genotyping, the discriminatory power of the panel 1 (MLVA-8) markers was less as compared to the combined panels 1+2A (MLVA-11) and 1+2A+2B (MLVA-16). All the Panel 2A and 2B primers have different HGDI values in different studies [27,28,57,58]. In our study, Bruce04 (panel 2B marker detecting 11 alleles) exhibited the maximum HGDI value, followed by Bruce16 (panel 2B marker detecting 9 alleles) and Bruce19 (panel 2A marker detecting 7 alleles) ( Table 1). The MLVA-16 analysis yielded the maximum number of MLVA genotypes. Similar results have been reported by other investigators [28,40,57,58]. A comparison with strains from other regions revealed considerable variation in the VNTRs associated with the same alleles [27,29,41]. Even the strains from same clade exhibit allele differences, except for Bruce 45, which represents a single allele in majority of cases [27,29,41].
In order to place our strains in a global perspective, we conducted the MST analysis to establish the genetic relatedness of the genotypes obtained in the current study with MLVA genotypes found worldwide. Based on MLVA-8 analysis, the most common genotypes (45 and 64), identified in the present study, were found in the East Mediterranean clade. Five out of seven novel genotypes also fell into the same clade. One novel genotype each belonged to the African and the West Mediterranean clades. The geographic influence plays a big role in the genotype distribution as none of the other common genotypes of the world were found in Kuwait (Fig 5). Kilic and co-workers reported that the strains isolated in Turkey were mostly from the East Mediterranean region and not from other parts of the world [29]. Interestingly, the genotype 64 was only at a distance of 1.0 with the genotype 43, which is commonly found in Turkey [29] and UAE [41]. Genotype 45 was also found previously in Turkey [39] and China [57]. Hence we assume that the Brucella have probably entered in Kuwait through animals and livestock imported from nearby regions. The strains isolated from Turkey [29], UAE [41], Lebanon [28] and China [57] also had genotypes belonging to the East Mediterranean region. The MLVA genotypes of Kuwait form a distinct branch in the East Mediterranean region with two separate groups. The process of evolution may have resulted in the formation of new genotypes. Further investigations on the novel genotype associated with the East Mediterranean region should be done. Owing to the further differentiation on the basis of panel 1+2a and Panel 1+2a+2b, the Kuwaiti arm of the East Mediterranean region forms a branched tree suggesting the presence of diverse genotypes of B. melitensis in the region.
In conclusion, our work demonstrates that the molecular techniques are fast and accurate tools for identifying and discriminating the strains of Brucella in Kuwait. The region is dominated by the pathogenic species of B. melitensis. 16S rRNA gene sequencing using MicroSeq 1 500 kit and real-time PCR can provide rapid confirmatory identification of Brucella isolates up to species level. The ERIC-PCR has a higher discriminatory power and a potent tool for intraspecies diversification. The technique suggests the presence of 75 ERIC genotypes in Kuwait. However, the ERIC-PCR is quite limited due to non-availability of comparative data set from other studies. The MLVA-16 genotyping scheme is capable of identifying the isolates up to species and genotype levels as well as trace back the origin of strains in a particular region. The strains in Kuwait have their origin from the East Mediterranean region and are in close resemblance with UAE strains.