Genome phylogenetic analysis of Brucella melitensis in Northwest China

Abstract

Brucellosis poses a severe threat to public health in Northwest China; however, the genome phylogeny and transmission pattern of Brucella melitensis from sheep and yaks in this region remain unclear. In this study, bacteriology, conventional biototyping, and whole-genome single-nucleotide polymorphism (WGS-SNP) were applied to depict the phylogenetic profiles of strains from Northwest China. A total of 46 Brucella strains were identified as B. melitensis bv. 3, which was isolated from at least three animal (livestock and wildlife) hosts, implying that B. melitensis infection is prevalent in the Northwest and suggesting that host diversity provides an optimal niche for the spread and maintenance of B. melitensis in this region. WGS-SNP analysis divided the 46 B. melitensis strains into four clades (C-I–IV) that harbored eight SNP genotypes (STs), implying that at least four lineages are prevalent in the Northwest. Global WGS-SNP phylogenetic analysis of strains revealed that all Northwest strains belong to genotype II. Strains from different clades presented high genetic similarity with strains previously collected from the Northwest. This study provides robust evidence supporting the notion that multiple similar B. melitensis lineages are persistently prevalent in human populations and animals in the Northwest. The economic development of animal husbandry has accelerated the cross-regional flow of livestock and livestock products, driving the spread and reach of the disease. Therefore, tailoring a targeted control strategy is necessary to counter the current serious epidemic trend.

Introduction

Brucellosis is a common zoonotic disease caused by Brucella spp. that severely threatens public health and leads to considerable economic losses in animals worldwide [1]. The burden that the disease places specifically on low-income countries has led the World Health Organization (WHO) to classify it as one of seven neglected zoonotic diseases [2]. The disease affects mainly livestock and wildlife, as well as food safety and international trade [3]. Animals contract the infection by ingesting contaminated feed and water or by contacting aborted fetuses or fetal membranes that are subsequently discharged from the uterus [4]. The disease can be transmitted from animal reservoirs to humans through direct or indirect contact with infected animals and contaminated material or through the ingestion of raw milk or unpasteurized dairy products [5]. An evidence-based conservative estimate of the annual global incidence is 2.1 million; Africa and Asia sustain most of the global risk and cases [6]. The incidence of human brucellosis has significantly increased in 31 mainland provinces in China, with a significant increase of 8.20% from 2004 to 2021 [7].

In China, animal brucellosis has experienced a persistent increase in the frequency of outbreaks and the number of reported cases from 2006 to 2021, with over 98% of reported cases occurring in sheep and cattle [8]. The annual incidence of brucellosis in humans is significantly positively correlated with the sheep population [9]. The pooled prevalence of brucellosis in ovine and caprine flocks in China increased from 1.00% from 2000 to 2009 to 3.2% from 2010 to 2018, and sheep populations are overwhelmingly predominant in major pasturing areas in mainland China, including Inner Mongolia, Xinjiang, Gansu, and Qinghai Provinces, which are categorized as severe epidemic regions of Type I brucellosis in China [10]. In these regions, B. melitensis is the dominant species in terms of prevalence and causes more than 84.5% of human brucellosis cases in China [11, 12]. The migration of humans and animals and the trade of animal products has created a challenge for disease spread among areas [13]. Many common and unique factors, including extensive husbandry, budgetary limitations, misdiagnosis, and other conditions, play a role in the long-term endemicity of brucellosis in these locations [14]. Although some animal vaccines, such as the B. abortus S19 and B. melitensis Rev. 1 vaccines, are the cornerstones of control programs in cattle and small ruminants, respectively, in the absence of a human brucellosis vaccine, the prevention of human brucellosis depends on the control of the disease in animals [15].

A retrospective spatiotemporal analysis of human brucellosis from 2019 to 2023 revealed that the largest cluster occurred in Northwest China [16]. MLVA genotyping of 46 B. melitensis strains from Northwest China revealed that strains from different regions may be epidemiologically linked [17]. B. melitensis bv. 3 was the dominant species in both humans and animals in Xinjiang and wgSNP results indicated that three main complexes were involved in the B. melitensis epidemic [18]. In Qinghai, WGS-SNP analysis of 54 B. melitensis bv. 3 strains revealed that cross-regional transmission events (i.e., between counties) were caused by common sources of infection [19]. However, the transmission pattern and genetic diversity of B. melitensis in Northwest China remain unclear. Single-nucleotide polymorphisms (SNPs) identified via whole-genome sequencing have allowed the characterization of the phylogenetic relationships of strains from different regions [20]. Therefore, this study aimed to investigate the transmission patterns and molecular relationships among these isolates using WGS-SNP phylogenetic analysis and to provide useful clues for devising tailored regional control measures.

Methods

Sample collection and isolation of strains

A total of 480 spleens from the aborted fetuses of multiple animal hosts (360 from sheep/goats, 108 from yaks, 9 from deer, 2 from cows, and 1 from a dog) were collected from the Northwest region of China, including the Gansu (n = 93), Inner Mongolia (n = 142), Xinjiang (n = 47), and Qinghai (n = 198) Provinces. Brucella was isolated as described in a previous study [17]. Briefly, 200 mg of each fresh and sterile spleen tissue sample was completely ground using a tissue homogenizer and then suspended in 500 µL of sterile phosphate-buffered saline (0.01 mol/L, pH 7.2). The prepared samples were subsequently uniformly coated on Brucella serum dextrose agar (containing Brucella Medium Base and Brucella Selective Supplement (OXOID, England) with 5–10% horse serum) in a biosafety cabinet according to sterile techniques [21]. The plates were subsequently inoculated with sample materials and incubated with or without 5–10% carbon dioxide at 37 °C. These plates were observed on alternate days for at least two weeks until transparent colonies appeared. A single pure suspected colony was selected for further species/biovar identification. Plates without transparent colonies after two weeks of incubation were regarded as negative and safely disposed of after autoclaving.

DNA preparation and species/biovars determination of strains

Total genomic DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen, Germany) according to the manufacturer’s instructions, and the prepared DNA samples were stored at –20 °C. Furthermore, AMOS–PCR was used to determine the species as previously described [22]. Briefly, PCR was performed in a 50 µL volume system containing 25 µL of Premix Taq, 1 µL of DNA template, and 0.2 µL of primers, with ddH₂O added to a final volume of 50 µL. The PCR amplification parameters were as follows: 96 °C for 5 min, 95 °C for 1 min, 60 °C for 2 min, 72 °C for 2 min, and 72 °C for 5 min. After this, 5 µL of each PCR product was separated via electrophoresis on a 1.5% agarose gel (w/v). DNA from the B. melitensis 16 M, B. abortus 544, and B. suis 1330 reference strains was used as a positive control, and the corresponding amplification band sizes were 731, 498, and 285 bp, respectively. Subsequently, agglutination with monospecific A and M antisera tests was used for biovar determination according to the standard test approach [23].

Genome sequencing

The genomic DNA was randomly sheared into short fragments, after which the obtained fragments were end repaired, A-tailed, and further ligated with an Illumina adapter. The fragments with adapters were size-selected, PCR amplified, and purified. The library was checked with a Qubit and real-time PCR for quantification and a bioanalyzer for size distribution detection. The quantified libraries were pooled and sequenced on Illumina platforms according to the effective library concentration and data amount needed. The original data were filtered to obtain valid data (clean data) to ensure the accuracy and reliability of the subsequent information analysis results. The data were assembled with SOAPdenovo [24], SPAdes [25], and Abyss [26] software. The assembly results of the three software programs were integrated with CISA software [27], and the assembly result with the least number of scaffolds was selected. Gap-closer software [28] was used to fill the gaps in the preliminary assembly results. Same-lane contamination was removed by filtering the reads with low sequencing depth (less than 0.35 of the average depth) to obtain the final assembly results.

Phylogenetic analysis of B. melitensis at the Northwest scale

Brucella melitensis 16 M (GCA_000007125.1) was used as a reference genome, core-genome SNPs (cgSNPs) were identified via snippy v 4.0, and snp-dists (v 0.6.3) was used to calculate the SNP distance between each pair of genomes. Maximum likelihood trees of 46 B. melitensis strains were generated via IQ-tree 2 using the maximum likelihood phylogenies algorithm with 1000 bootstrap replicates [29]. Furthermore, a genetic relationship compared with 110 strains was conducted with above-mentioned algorithm methods, of which 64 (47 previously collected from the Northwest, and the remaining 17 from five (I–V) different genetic lineages [30]) were downloaded from the GenBank database (Table S1). The tree was visualized and edited using iTOL (Interactive Tree Of Life) v6.5.7 [31].

Results

Species/biovars and distribution of 46 Brucella strains

AMOS-PCR amplification revealed a special band of 731 bp in size in all 46 strains, and the strains were B. melitensis. Furthermore, the results of a positive serum agglutination test with anti-M and anti-A sera revealed that the 46 strains were all B. melitensis bv. 3 (Table 1). The strains were from four Northwestern provinces: 21 from Inner Mongolia, 13 from Xinjiang, 11 from Gansu, and 1 from Qinghai (Table 1), of which 42 were from sheep, 2 from deer, 1 from humans, and 1 from yak. These data indicate that B. melitensis bv. 3 is widely prevalent in at least four hosts in the Northwestern region.

Table 1 Distribution of location, hosts, isolation time, species/biovars, STs, and number of strains involved in this study

Epidemiological correlation analysis of 46 Brucella strains

Based on WGS-SNP analysis, 46 B. melitensis strains were divided into four clades (C-I–IV) that harbored eight SNP genotypes (STs) (Fig. 1). C-I and C-II each consisted of one ST, C-III was composed of two STs (ST3 - 4), and C-IV harbored four STs (5–8) (Fig. 2). In addition to C-I, strains from each clade were isolated from at least two provinces: C-II included strains from Inner Mongolia and Gansu; C-III included strains from Inner Mongolia, Gansu, and Xinjiang; and C-IV included strains from four regions, including Inner Mongolia, Gansu, Xinjiang, and Qinghai (Fig. 2). Two strains from deer presented the same SNP genotype as one strain from sheep from Xinjiang; one strain from a human from Gansu presented high genetic similarity with two strains from sheep from Inner Mongolia; similarly, one strain from yak in Qinghai was closely related to one strain from sheep from Inner Mongolia. In addition, strains from two STs (6 and 8 (partly)) contained strains from a single region (Inner Mongolia and Xinjiang) (Fig. 1). These data implied that there was a co-epidemic pattern of animal brucellosis across regions and local transmission in the four Northwestern provinces.

Fig. 1

Maximum-likelihood tree generated from the cgSNP matrix of 46 B. melitensis samples at the county scale. Note: The phylogenetic trees of 46 B. melitensis strains were generated via TreeBeST via the maximum likelihood phylogenies (PHYML) algorithm with 1000 bootstrap replicates. Note: A branch marked by color indicates the location of strains isolated; the red branch indicates that strains are from Gansu; the blue branch indicates that strains are from Qinghai; the green branch indicates that strains are from Inner Mongolia; and the gray branch indicates that strains are from Xinjiang. The host spectrum is indicated in the right color column; the red column indicates strains from sheep, the blue column indicates strains from humans, the green column indicates strains from deer, and the yellow column indicates strains from yaks. ST1 - 8 indicates that strains are associated with each SNP genotype

Fig. 2

Phylogenetic analysis based on the maximum-likelihood tree of 110 B. melitensis strains at the global level. Note: The strains from GTs I–III refer to previously described methods [30] and are marked in red. The phylogenetic trees of 110 B. melitensis strains were constructed based on the maximum likelihood phylogenies algorithm with 1000 bootstrap replicates. Note: Among the 110 strains, 15 strains from GenBank display the five SNP genotypes (I–V, marked in red) of B. melitensis, and the remaining 95 strains were all isolated from Northwest China from GenBank. Strains from this study are marked with bold italics, and branches marked by color indicate the locations of the strains isolated; the green branch indicates that the strains are from Inner Mongolia; the red branch indicates that the strains are from Gansu; the blue branch indicates that the strains are from Qinghai; and the purple branch indicates that the strains are from Xinjiang. The host spectrum is marked in an inward circle, with humans in blue, gray in goats, dark green in cattle, red in sheep, light purple in cows, green in deer, and yellow in yaks. The years in which the strains were isolated are marked with external circles; the strains from 1957 to 1990 are marked with light yellow, the strains from 2007 to 2015 are marked with light blue, and the strains from 2016 to 2022 are marked with dark green

Phylogenetic analysis of 110 B. melitensis from multiple hosts in the Northwest

A global comparison analysis revealed that all strains in this study belonged to genotype II, and strains from Northwest China were closely related to strains from genotypes II f–i (Fig. 2). Strains from ST1 belonged to genotype II f; strains from Xinjiang and Gansu formed a subclade with strains from Geo125210 (genotype II i) isolated from Georgia; strains from ST2 presented high homogeneity with strains from Xinjiang; strains from ST3 and ST4 clustered with strains from Inner Mongolia and Xinjiang; strains from ST5 and ST6 clustered with strains from Inner Mongolia; ST7 grouped with strains from Inner Mongolia and Xinjiang; strains from ST8 presented high homogeneity with strains from Inner Mongolia. These data revealed that strains from the Northwest had no apparent geographical boundaries; furthermore, strains from different hosts and isolation times shared the same SNP genotype (Fig. 2). These data implied that similar B. melitensis lineages, including humans, sheep/goats, deer, and yaks, were persistently spread and expanded in Northwest China in multiple hosts.

Discussion

Although brucellosis has been successfully managed in most industrialized countries, it significantly burdens on goat and human health in the Mediterranean region, the Middle East, Central Asia, and Southeast Asia, including China [32]. In this study, a genome phylogenetic analysis of strains from Northwest China was performed to better understand the epidemiological profile. The 46 B. melitensis strains were isolated from abortus fetuses from three animal hosts, and all the strains were identified as B. melitensis bv. 3 These data suggest that B. melitensis was the dominant species in multiple hosts in Northwest China. From 2010 to 2016, a total of 110 B. melitensis strains were obtained from tissue and milk samples collected from cattle and small ruminants at multiple locations in 7 northern provinces and 50 counties [33]. This conclusion is supported by 20 osteoarthritis-associated B. melitensis strains isolated from inpatients in Inner Mongolia from 2013 to 2017 [34]. Twelve Brucella strains, which are recognized as B. melitensis, were identified in 62 patients in Xinjiang [35]. There are at least three Brucella (B. melitensis, B. abortus, and B. suis) species in Qinghai, of which B. melitensis is the predominant species in the area examined [36]. Sichuan borders Northwest China; a total of 101 Brucella strains were isolated from 16 cities (autonomous prefectures) from 2014 to 2021, and all the strains were identified as Brucella melitensis bv. 3 [37]. In Northwest China, the populations of sheep and goats account for an overwhelming majority of the populations in other regions, and the prevalence rates of brucellosis in ovine and caprine flocks significantly increased from 1.00% between 2000 and 2009 to 3.20% in 2010 and 2018 [10]. In Inner Mongolia, the largest number of human brucellosis cases occurs in mainland China, which can be attributed to the large number of sheep kept there; at least 90% of human brucellosis cases are caused by sheep [38]. In 2021, the number of counties with an incidence of human brucellosis greater than 10.00/100,000 was 517; these counties were located in Inner Mongolia (95), Shanxi (79), Xinjiang (67), Heilongjiang (56), Hebei (38), and Gansu (37) [39]. These data indicate that B. melitensis is a regional public concern, and more resources need to be allocated to Northwest China to strengthen the prevention and control of animal brucellosis to prevent its spread.

A global phylogenetic analysis revealed that all the strains in this study belonged to genotype II, and the strains presented high genetic similarity with strains from Georgia, Pakistan, and Kuwait. The strains of genotype II have the widest geographical distribution, and their habitat is most of the territory of the Eurasian continent, from Cyprus in the west to the islands of Southeast Asia in the east, suggesting a possible single introduction of Mediterranean origin [30]. The B. melitensis East Mediterranean lineage is predominant across the continent, with only a small number of samples from the Africa and Americas lineages [40]. The strains from ST1 had high genetic similarity with B. melitensis BG2 (S27), which was isolated from Pakistan and occupied the ancestral node of group f (China clade), suggesting that the B. melitensis strain currently circulating in China may have a Pakistani origin [41]. Frequent animal trade and exchange have provided opportunities for B. melitensis prevalence in these Asian countries [42].

Furthermore, WGS-SNP analysis revealed the absence of clear transboundary differentiation according to strain territorial affiliation, implying that similar or identical B. melitensis lineages were persistently spread and expanded in Northwest China. In addition, these B. melitensis lineages are shared by strains from multiple different hosts, including humans, livestock (sheep/goat), and wildlife animals (deer and yak); an extensive host spectrum provides an optimal niche for the spread and maintenance of B. melitensis in Northwest China. On the Qinghai‒Tibet Plateau, B. melitensis strains isolated from humans and marmots exhibit high genetic similarity, implying the possible existence of direct or indirect contact between sheep (and goats) and wildlife (marmots) [19]. Whole-genome sequencing of strains from outbreaks of brucellosis in elk in the Greater Yellowstone Ecosystem demonstrated that free-ranging elk are currently a self-sustaining brucellosis reservoir and a source of livestock infections and that control measures in bison are unlikely to affect the dynamics of unrelated strains circulating in nearby elk populations [43]. Multiple host infections caused by B. melitensis strains pose a great public health challenge for the surveillance and control of this disease. In addition, frequent animal husbandry and economic trade have accelerated the cross-regional flow of livestock and livestock products, and the inadequate implementation of quarantine measures for transregional livestock transport; these factors have co-driven the spread and expansion of the disease. Brucellosis is associated with large-scale farming and trading of sheep and goats; the nonregulated and frequent animal trade greatly contributes to the dissemination of brucellosis [11]. Infected sheep/goats and their products before import from north to south may be the main cause of the increasing incidence of human brucellosis in southern China [44]. Owing to population migration and trade transport, brucellosis continues to prevail in Asia and Africa and is now affecting regions previously considered under control, such as the United States [45]. The movement and illegal trade of infected sheep and their products require stringent control to reduce the risk of human infection.

Conclusion

In the present study, B. melitensis bv. 3 is widely prevalent in different animal hosts in the Northwestern region, which poses a serious threat to public health and implies that the surveillance and control of brucellosis in this region should focus on B. melitensis strains. Further WGS-SNP analyses of strains from different regions and hosts revealed high genetic similarity, and these data provide robust evidence to support the notion that similar B. melitensis lineages are persistently circulating in multiple hosts in Northwest China. These results provide vital information for the regional and national surveillance and control of brucellosis. A molecular analysis of strains on a nationwide scale will provide comprehensive insight, which will allow us to better understand the transmission dynamics and genetic variant profile of B. melitensis.