Genomic Characterization of Clinical Brucella melitensis Isolates from Tunisia: Integration into the Global Population Structure

1. Introduction

Human brucellosis, a disease caused by the Brucella species of bacteria, is one of the most prevalent zoonotic diseases globally, with approximately 500,000 new cases occurring annually [2]. The epidemiology of this disease has drastically changed over the past decade because of various sanitary, socioeconomic and political reasons [1]. Despite the availability of control measures such as animal vaccination, surveillance of abortions and pasteurization of milk products, brucellosis remains endemic in many countries of the Mediterranean region, the Middle East, Southern Europe, Central Asia, the Indian subcontinent, and Africa [1,3]. In 2017, the overall incidence in Tunisia was 9.8 per 100,000 populations with 621 new human cases in the mentioned year [4].

Humans represent an occasional host infected via direct contact to infected sheep or goats or indirectly through the consumption of unpasteurized milk, raw meat, and dairy products contaminated with the causative bacteria [5]. Therefore, brucellosis represents a significant public health concern as well as an economic burden as compromised food security requires counteractive measures (quarantine or slaughter) with critical impact on livestock productivity [1]. Human-to-human transition is rare, it can occur transplacentally, via breastfeeding, and rarely through sexual intercourse, organ transplantation and blood [3]. The Brucella genus comprises one of the most common pathogens causing laboratory-acquired infections (LAIs) which are becoming an increasingly important biosafety issue [6]. Brucella species are also considered potential bioterrorism agents [7]. If acquired and adequately disseminated, Brucella could cause a serious public health challenge in terms of the ability to limit casualties and control other consequences of such an attack [8].

The genus Brucella comprises more than 12 species with a percentage of similarity close to 100% through DNA-DNA hybridization [9,10]. Of these, B. melitensis is the predominant species responsible for the majority of human cases [11,12,13]. Despite their close genetic relationship, Brucella species exhibit a wide diversity of phenotypic characteristics, host preferences, and pathogenicity thus requiring robust universal diagnostic tools. A high discrepancy between the reported rate and the actual incidence has been reported largely and it is due to misdiagnosis and underdiagnosis, especially in endemic areas [14]. Currently, the diagnosis of human brucellosis in Tunisia is based on a combination of clinical examination, serological tests, cultivation and bacterial isolation from biological samples. These methods are generally time-consuming and required laboratory-qualified personnel, as well as a biosafety level 3 facility [15,16]. If laboratories are equipped for it, real-time PCR based assays are the method of choice for the rapid and reliable diagnosis of brucellosis. The assays currently in use focus on identifying different species of the genus Brucella based on its highly conserved genetic loci (e.g., BCSP 31, IS6501/711 or 16S rRNA genes).

Core genome multilocus sequence typing (cgMLST) is widely used to assess genetic relatedness among bacterial isolates to support outbreak investigation [17,18] or tracing back transmission chains [19]. Whole Genome Sequencing (WGS) based phylogeny stands as the most powerful tool to discriminate Brucella strain identity, origin and their phylogenetic relationship [20]. This approach facilitates a comprehensive understanding of the global epidemic.

However, the lack of epidemiological data from Tunisia and its neighboring countries represents a significant obstacle for the assessment of the clinical impact of this region. The objective of this study was twofold: Firstly, to examine the classification, within the context of both the local epidemiological landscape and the global epidemiological picture, of strains identified in Tunisia through the use of whole genome sequencing. And secondly to describe virulence patterns of these clinical Brucella melitensis strains isolated in Tunisia.

2. Materials and Methods

2.1. Isolation of Bacteria

A total of 36 strains of suspected Brucella spp. were obtained from 36 blood cultures collected between January 2016 and December 2020 from patients coming from different regions of Tunisia. They were suspected of acute brucellosis and hospitalized in the infectious diseases departments of two university hospitals in Tunisia: Hôpital Charles Nicolle (the biggest and the oldest hospital in Tunisia) and Hôpital La Rabta. The latter being the reference center of infectious diseases responsible for handling patients in need coming from different northern regions. Blood samples were placed in an automatic incubator (Bactalert 3D, BioMérieux, Marnes La Coquette, France). Positive blood cultures were plated on Columbia blood agar (BioMérieux) containing 5% sheep blood and incubated for 48 h at 36 ± 1 °C with 5% CO₂. All strains were studied for colony morphology and culture features and were characterized by morphological (e.g., microscopic stain) and biochemical tests such as oxidase, catalase, urease and Api 20E gallery identification (BioMérieux, Marnes La Coquette, France). Strains identified as Brucella spp. were stored at -80°C in Brain Heart Infusion (BHI) with 15% glycerol.

2.2. Identification of Bacterial Strains Through Real-Time PCR

Total DNA was extracted from single colonies of subcultures using a DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol for Gram-negative bacteria. For the multiplex PCR (table 1), two Brucella gene markers were targeted (IS711 and MazG), as well as an internal positive control (IPC). The IPC is a chemically synthetized sequence (Koma2) cloned into the pPCR-Script Vector (Stratagene, LA Jolla CA). It is added to the samples before extraction at a concentration of 250 copies per sample. Detection of its presence at an average Ct value of 35 ± 3 within the multiplex PCR verifies the success of the extraction and the functionality of the PCR set up. The localization of the primers and probes for Brucella detection is based on the sequence of the MazGgene (GenBank U78089) and the IS711 gene (GenBank M94960). PCR amplification was performed in a final volume of 25 µL. The reaction mixtures consisted of 5 µL of DNA template, 6.25 µL TaqMan Environmental Master Mix 2.0 (Life Technologies, Carlsbad, USA), 0.75 µL of each primer (10 µM) and 0.25 µL of probe (10 µM) and 8.5 µL of free nuclease water (Fluka, Lincolnshire, IL, USA).

For amplification, an ABI 7500 system (Applied biosystems, Waltham, MA, USA) was used following a standard amplification program (table S1). DNA from the reference strain of Brucella melitensis bv. 1 str. 16M was used as a positive control. The concentration of this DNA was adjusted to 0.2 pg µL^-1 by dilution with AE buffer (Qiagen, Hilden, Germany) so that a Ct value between 31-35 (depending on the target marker) was obtained by real-time PCR.

2.3. Whole Genome Sequencing and Annotations

Total genomic DNA was quantified with the Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA), and library preparation was performed using the Nextera XT library preparation kit (Illumina, San Diego, CA, USA). The libraries were sequenced using the Illumina NextSeq 500 platform, producing 150-bp paired-end reads. After demultiplexing and removal of adapters, reads were trimmed and filtered with fastpv0.22.0 using default options [21]. Quality control of raw and filtered reads was conducted with FastQCv0.11.9 [22]. The de novo whole-genome assemblies were performed with SPAdesv3.15.5 using the option –careful [23]. Contigs with a length under 1000 bp were removed with bbmap v39.01. The species Brucella melitensis and the percentages of completeness and contamination were inspected with CheckM v1.2.2 [24]. Virulence and antimicrobial resistance (AMR) genes in all Tunisian strains and two reference strains were detected using ABRicate v1.0.1with the databases VFDB, and MEGARes and AMR FinderPlus3.11.4 with the NCBI database and the option –plus.

2.4. Core Genome Multi-Locus Sequence Typing (cgMLST)

A cgMLST scheme for B. melitensis with 2704 targets in the core genome and 360 targets in the accessory genome was used [25]. Basically, a genome-wide gene-by-gene comparison using the cgMLST Target Definer (version 1.4) function of the SeqSphere software, v5.0.90 (Ridom GmbH, Münster, Germany), was performed. Default parameters and the following filters were selected to exclude certain genes of the B. melitensis bv. 1 strain 16M reference genome (NC_003317.1 and NC_003318.1) from the cgMLST scheme: a minimum length filter that discards all genes shorter than 50 bp, a start codon filter that discards all genes that contain no start codon at the beginning of the gene, a stop codon filter that discards all genes that contain no stop codon or more than one stop codon or if the stop codon is not at the end of the gene, a homologous gene filter that discards all genes with fragments that occur in multiple copies within a genome (with identity of 90% and more than 100-bp overlap), and a gene overlap filter that discards the shorter gene from the cgMLST scheme if the two genes affected overlap by 4 bp. The remaining genes were then used in a pairwise comparison using BLAST (version 2.2.12), with the query chromosomes of one representative for each of the other two B. melitensis biovars (B. melitensis bv. 2 strain 63/9 [NZ_CP007788.1 and NZ_CP007789.1] and B. melitensis bv. 3 strain Ether [NZ_CP007761.1 and NZ_CP007760.1]). Using all genes of the reference genome that were common in all query genomes, with a sequence identity of 90% and 100% overlap, and with the genome filters start codon filter, stop codon filter, and stop codon percentage filter turned on, the final cgMLST scheme was formed. Therefore, all genes having no start or stop codon in one of the query genomes, as well as genes that had internal stop codons in more than 20% of the query genomes, were discarded. An additional cgMLST was performed using the pipeline and database of PubMLST for Brucella species [26]. The cgMLST scheme from Tunisian strains was obtained using the tool Genome Compartor in PubMLST.

2.5. Clustering Analyses, Cluster Definition, and Data Availability

The cgMLST profiles were assigned using B. melitensis task template in SeqSphere. Minimum spanning trees (MSTs) were created by a pairwise comparison of cgMLST target genes using default software settings. Missing values were ignored in the calculation of the distance between pairs of sample profiles. The links between the MST nodes represented the distance between the genotypes. The cluster cutoff value was defined as the maximum pairwise distance found between epidemiologically related isolates. In the case of PubMLST, MSTs were created with GrapeTree using default parameters [27].

All generated data was submitted to the National Center for Biotechnology Information (NCBI) under the Bioproject IDPRJNA1053259. The 36 assemblies were assigned an accession number from SAMN38851271 to SAMN38851306 (see table S2 for further details).

2.6. Acquisition of Metadata

Demographic and clinical data were obtained from human brucellosis cases infected by the isolated strains. Physicians from different department were contacted to complete an information sheet for each patient including: origin, age, gender and symptoms, (see in table S2).

3. Results

3.1. Collection of Isolates and Molecular Identification

Respective demographic and clinical data on the 36 Brucellosis cases are summarized in Supplementary Table S2. The real-time PCR assays performed in parallel at the Robert Koch Institute (RKI) in Berlin (Germany) and at Charles Nicole Hospital in Tunis (Tunisia) yielded identical results: all isolates were identified as B. melitensis.

3.2. Whole-Genome Analysis and Core Genome Multi Locus Sequence Typing (cgMLST)

All cultivable isolates were subjected to whole genome sequencing. All assemblies yielded an average coverage of ~130x and corresponded to B. melitensis with over 99% completeness and less than 1.3 % contamination. The neighbor-joining tree shows that, in congruence with PCR-based species identification, all Tunisian isolates cluster together with B. melitensis biovar 3 (Fig. 1A).

Both tested schemes used to determine genetic relatedness between the isolates were comparable and revealed identical groupings with only few allele differences between the two assays. Within the global population structure, Tunisian isolates form two distinct clusters with a distance of 140 discriminating alleles. Within the proposed clusters, strains differed in up to 6 alleles, with 14 isolates from proposed Tunisia cluster 1 showing not a single discriminatory allele. In addition, the MST shows that the Tunisian clusters are most closely related to countries in geographical proximity, with 150 discriminating alleles to Egypt, 245 to Algeria, and 266 to Morocco (Fig. 1B and 2). Of the 36 analyzed strains 5 did not belong to any cluster, with one of them being very closely related to proposed Tunisia cluster 1. One of these isolates appears to be most closely related to single isolates from Morocco and Algeria. Interestingly, with less than 30 discriminating alleles, isolates from an outbreak in Austria cluster very close to the proposed Tunisia Cluster 2 and 3 other strains from this sample set which could not be allocated to a cluster.

3.3. Identification of AMR and Virulence Genes from WGS Data

All Tunisian strains presented the same virulence and AMR genes than the reference strains of the biovars 1 and 3 (Table 2). For example, type IV secretion system genes (e.g, virB, ricAand bsp), immune modulation (lpx) and adhesins (big) were in every strain detected.

4. Discussion

This study is one of the few Tunisian studies investigating human Brucella strains. Brucellosis is a zoonosis of worldwide distribution that is endemic in the Mediterranean basin, especially in the Northern African countries [28]. Despite the important animal reservoir for Brucella species in Tunisia, the epidemiological situation of human brucellosis was not known in Tunisia until the 1980s, particularly because of the misdiagnosis of the disease [29]. Since the beginning of the 1990s several outbreaks have occurred, such as the Gafsa (southern Tunisia) outbreak in 1991, with more than 400 human cases posing a serious public health threat. Hence eradication and control measures in both human and animal brucellosis were conducted, in Tunisia, the overall incidence for human brucellosis ranged between 2.9 and 4.5 per 100.000 inhabitants in 2015 and 2023, respectively. Which can explain the significant number of cases recorded in 2016 and 2017 [30]. Despite the implemented programs, Brucellosis is still occurring in Tunisia, with a variable annual incidence. As a chronic disease, it is of particular importance for public health, with indirect economic impacts such as decreasing the work force and reducing livestock reproduction and their associated products [31,32].

In this study, Tunisian human Brucella strains were clearly identified and characterized. A good and precise diagnosis with an exact identification of the bacterial species is essential to knowing the disease incidence in a region, to planning epidemiological studies, and to controlling eradication programs. Our results support the previous finding that human brucellosis can be a common public health problem in many governments in Tunisia [33]. A comprehensive picture of the epidemiology and genetic diversity of B. melitensis in Tunisia was provided.

Tunisian B. melitensis strains presented an elevated genome similarity and shared almost all virulence and AMR-related genes, even the same genes as the considered reference strains. Brucella species are well-known to present almost identical genomes in comparison with other bacterial species [25].

One of the strengths of this study lays in the analysis of a database data on Human B. melitensis strains; this approach allows more comprehensive understanding of the local situation regarding human brucellosis in Tunisia. However, a limitation of this study is the number of strains, emphasizing the need for larger datasets for more accurate future research.

Our phylogenetic analysis of B. melitensis followed similar phylogeographical patterns, and two main phylogenetic clusters were highlighted, with Tunisian Cluster 1 being more clearly defined (25 isolates) by the number of Isolates than Tunisian cluster 2 (6 isolates). Isolates from both populations were found essentially to circulate in the north and the north west of Tunisia. The overall distribution in Tunisia cannot be assessed due to the bias of the catchment area of the two hospitals involved in the sampling. Of note, in 2016 and 2018, two important outbreaks of human brucellosis were observed in the south of Tunisia [30]. The two main branches closely related to Cluster 1 were linked to strains from Italy (the closest) and Egypt. The history of brucellosis spread in the Mediterranean basin could be explained by the strong trade connections between Italy and North Africa, lasted from the time of the Roman Empire to the Middle Ages, and continued in the last century, during the period of modern colonialism [25]. Commercial activities between both regions could explain the circulating Brucella species between African countries and the Northern Mediterranean region [29]. This study evidenced that the Austrian isolates were clustered together with Tunisian strains, indicating an association between both countries. Austrian isolates belonged to a branch including Tunisian Cluster 2. The Tunisian cluster 2 and were grouped together with Italian clusters, suggesting the hypothesis that both Italian and Austrian isolates were imported from Tunisia. This is especially supported by figure 2, where a smaller Italian branch is calculated to originate from the strains identified in this study. A relationship between these countries was previously described in an Austrian study, which showed that B. melitensis isolates collected between 2005 and 2015 were closely related to a branch including strains from Egypt and Italy [30]. The studies characterizing the population of B. melitensis using WGS remain limited to some Mediterranean countries of the Maghreb region (Lybia, Algeria, Morocco) [1], we included recently published WGS data from Tunisia to expand the phylogeny and provide a comprehensive analysis.

To date, a few strains have been isolated from Africa, and the lack of B. melitensis genomes from other European countries within the Mediterranean region prevents us from knowing the real history of brucellosis in the region [28].

Implementation of high-throughput WGS to identify the AMR and virulence-associated genes in Brucella isolates revealed no apparent difference in their distribution between B. abortus and B. melitensis strains or isolates from different hosts [34]. Investigations on antibiotic susceptibility and updating of breakpoints are required. Investigations of resistance and virulence mechanisms at the proteomic and transcriptomic levels have to be considered in future research.

In conclusion, this study has helped to fill one of the gaps regarding the phylogenetic nature of the strains circulating in Tunisia in relation to its neighbouring countries. This is an important step towards understanding the evolution of Brucella melitensis in this region. However, more data are needed from southern Tunisia and its neighbouring countries to complete the picture and allow a more confident interpretation of the phylogenetic relationships presented today.