Genetic diversity atlas of Brucella melitensis strains from Sichuan Province, China

Abstract

Human brucellosis is a re-emerging disease in Sichuan Province, China. In this study, bacteriology, conventional bio-typing, multi-locus sequence typing (MLST), and multiple locus variable-number tandem repeat analysis (MLVA) were applied to preliminarily characterize the strains in terms of genetic diversity and epidemiological links. A total of 101 Brucella strains were isolated from 16 cities (autonomous prefectures) from 2014 to 2021, and all of the strains were identified as Brucella melitensis bv. 3, suggesting that surveillance should focus on ruminants. MLST analysis identified four STs, namely, ST8 (n = 93), ST39 (n = 6), ST101 (n = 1), and ST118 (n = 1). The latter were new STs, indicating that strains displayed high population diversity. Six MLVA-8, namely, 42, 43, 45, 63, 83, and 114, and eight MLVA-11, namely, 111, 115, 116, 125, 180, 291, 298, and 342, genotypes were identified, demonstrating that all of the strains were from the Eastern Mediterranean lineage, and these strains exhibited a high genotype diversity. MLVA-16 analysis revealed that there was a co-existing transmission pattern, where sporadic cases and multiple outbreak events had a common origin. The dominant STs and MLVA genotypes of strains were epidemic in Northern, China, and 36 MLVA-16 genotypes were shared among strains (n = 51, 50.4%, 51/101) from Sichuan and strains from 22 other provinces. The findings imply that infected animals were introduced from outside the province. The surveillance and control of the disease have become public health challenges. Animal quarantines should be strengthened to prevent the spread of B. melitensis species among adjacent regions.

Introduction

Brucellosis, a globally widespread zoonosis, is caused by Brucella spp., a facultative intracellular pathogen. Brucellosis is a reproductive disease in domestic animals and a chronic debilitating disease in humans [1]. Humans, especially in low-income countries, can acquire infection by direct contact with infected animals or by the consumption of animal products. The disease can cause severe economic losses and public health problems [2]. Although B. abortus brucellosis has been eliminated or brought under control in developed countries [3], the disease remains one of the most common bacterial zoonoses in many low-income areas, including Africa and Asia. The highest incidence of human brucellosis is in countries in Asia [4].

In China, brucellosis is a re-emerging disease that expanded from northern pastureland provinces to the adjacent grassland and agricultural areas, then to the southern coastal and southwestern regions [5]. At present, human brucellosis has been reported from all 31 provinces or autonomous regions in mainland China. B. melitensis is the main pathogenic bacteria causing brucellosis epidemics in China [6]. Sichuan Province is located in the interior of Southwest China, bordering on Guizhou, Qinghai, Gansu, and Shaanxi provinces, and brucellosis is endemic in the region [7]. Sichuan Province was an endemic area of brucellosis from the 1950s to 1980s; the disease was brought under control in the 1990s and re-emerged in 2008 [8]. Since 2008, the reported number of cases has gradually increased. A total of 417 new cases were reported in 2023, and these were distributed among multiple cities. However, the level of genetic diversity and the epidemiological links of Brucella isolated from different regions and at different times need to be clarified. Therefore, bio-typing was used to determine the species/biovars of circulating Brucella strains, and the genotyping methods of multi-locus sequence typing analysis (MLST) [9] and multi-locus variable number tandem repeat analysis (MLVA) [10] were applied to investigate the genetic diversity and molecular links among the strains. The results provide valuable for the planning of a control strategy for brucellosis in Sichuan Province.

Materials and methods

Bacteriological isolation and bio-typing

A total of 101 Brucella strains were isolated from human blood samples obtained during the period from 2014 to 2021 using standard bacteriological methods [11]. The Rose Bengal plate test and serum agglutination test were used for brucellosis screening and diagnosis [12]. The procedures of blood collection, culture, incubation, and bio-typing were according to previously reported methods [13]. B. suis S2 vaccine strain was used as controls in the bio-typing. Among the 101 strains, there was one in 2014, seven in 2015, three in 2016, five in 2017, six in 2018, 14 in 2019, 28 in 2020, and 37 in 2021.

DNA extraction and genotyping of Brucella strains

All of the strains were inactivated at 100 °C for 10 min, and DNA was prepared using a Qiagen bacterial genome extraction kit (Qiagen, Heidelberg, Germany) according to the manufacturer’s protocol. The strains were verified to be B. melitensis through the detection of a 731 bp band specific to this species by AMOS-PCR [14]. Subsequently, based on the protocol of the study, the 101 strains were subjected to MLST and MLVA to quantify genetic diversity. Nine genomic loci were selected for MLST genotyping based on a previous study [9], namely, gap, aroA, glk, dnaK, gyrB, trpE, cobQ, omp25, and int-hyp. PCR amplification of the nine loci was performed as described previously [15]. Positive PCR products were evaluated by 1.0% agarose gel electrophoresis, and the products were sequenced. The sequences obtained from the purified PCR products were aligned using MEGA 6.0 software according to published allele MLST sequences (https://pubmlst.org/brucella/). Finally, the profiles of the nine loci were identified as specific sequence types (STs) based on the MLST online database (https://pubmlst.org/data/). MLVA genotyping referred to a previously described method [15]. In all, 16 loci were used to characterize the strains; these comprised three panels: panel 1 (MLVA-8), MLVA-11, and MLVA-16. Briefly, the PCR products were evaluated by 2% or 3% agarose gel electrophoresis, and the products were sequenced by capillary electrophoresis on an ABI Prism 3130 automated fluorescent capillary DNA sequencer (Applied Biosystems). Gene Mapper v. 4.0 (Applied Biosystems) was used to convert the fragment sizes of each locus to repeat unit numbers. A phylogenetic analysis (Table S1) was performed using the unweighted pair group method with arithmetic mean (UPGMA) in Bionumerics software v. 8.0 (Applied Maths, St-Martens-Latem, Belgium). According to a previous study [16], the locus diversity and discriminatory power of each method were analyzed based on the Hunter-Gaston diversity index (HGDI). To investigate the source of infections, the genetic relationships of 801 isolates (700 previously reported and 101 from this study) (Table S2) were evaluated, and the identical MLVA-16 genotypes were displayed using a minimum spanning tree (MST) constructed by PHYLOViZ 2.0 software [17]. A global genetic correlation analysis of 1296 B. melitensis strains (1195 B. melitensis strains from the international microbes genotyping MLVA database (http://microbesgenotyping.i2bc.paris-saclay.fr/databases) and the remaining 105 strains from the present study) (Table S3) was conducted using the minimum spanning tree (MST) based on Bionumerics software v. 8.0 to explore the genetic correlations among strains.

Results

Geographic and species/biovars distribution of Brucella strains

A total of 101 Brucella strains were isolated from samples taken between 2014 and 2021 (Table 1). Based on bio-typing approaches, all of the strains were identified as B. melitensis bv. 3. All of the strains possessed a specific 731 bp band amplified by AMOS-PCR. The B. melitensis strains were the predominant circulating species in this region. All of the strains were isolated from the blood of humans, were distributed in 16 cities, with the majority originating from Luzhou City (n = 19), followed by 17 from Neijiang City, 11 from Zigong City, and 10 from Chengdu City. The number of strains in the other 13 areas ranged from 1 to 9 (Table 1; Fig. 1).

Table 1 The distribution of genotypes in 101 B. melitensis strains

Fig. 1

The distribution profiles of strains, STs, and MLVA-11 genotypes

Note: The red dot (upper right figure) indicates the location of Sichuan in China

VNTR and allele diversity profiles

In MLST, the highest number of alleles was associated with the glk locus (n = 3), which also had the highest diversity, followed by trpE with two alleles, and the remaining seven loci with one allele ach (Table 2). The partition for MLVA was the highest for Bruce09 with 10 allele types. The allele numbers for Bruce04, Bruce16, and Bruce30 by MLVA were 8, 8, and 5, respectively, and the range of Simpson’s ID was from 0 to 0.808 for Bruce04, followed by 0.741 for Bruce16 and 0.671 for Bruce30. Moreover, three loci from panel 1 displayed low diversity; the allele numbers for Bruce12, Bruce42, and Bruce43 were 3, 2, and 3, respectively (Table 2).

Table 2 Locus diversity profiles and resolution of genotyping methods

Table 3 Alleles, number, and distribution profiles of STs considered in this study

MLST genotyping characteristics of 101 B. melitensis strains

Based on nine loci in the MLST assay, the 101 Brucella isolates were divided into four sequence types (STs): ST8 (n = 93; 92.00%), ST39 (n = 6; 5.94%), ST101 (n = 1; 0.99%), and ST118 (n = 1; 0.99%), of which ST101 and ST118 (n = 1) were newly isolated STs in this study; both were represented by one strain (Table 2). The ST8 strains were widely distributed among the 17 cities/states; the ST39 strain was distributed in six cities/states, while ST101 and ST118 only occurred in Leshan and Luzhou, respectively (Fig. 1). ST101 differed from ST8 in a single locus variant of the glk gene, while ST118 had seven loci that differed from ST8 (Table 3). The data indicated that ST8 was the dominant circulating Brucella species in Sichuan, and that the rise of new STs illustrates the need for active surveillance.

MLVA genotyping characteristics of 101 B. melitensis strains

In panel 1 (MLVA-8), all of the strains comprised six genotypes, namely, 42 (n = 72, 71.00%), 43 (n = 24), 45 (n = 1), 63 (n = 1), 83 (n = 1), and 114 (n = 2). Seven MLVA-11 genotypes were identified: 111 (n = 1), 115 (n = 1), 116 (n = 70, 69.00%) and 125 (n = 24), 180 (n = 1), 291 (n = 2), 298 (n = 1), and 342 (n = 1), with 116 being the dominant circulating genotype in the CI cluster and 125 dominating the CII cluster. The results suggest that all of the strains belonged to the East Mediterranean lineage. Based on MLVA-16 technology, the 101 strains were grouped into two clusters (CI and CII) and 13 sub-clusters (a ~ m) that harbored 74 MLVA-16 genotypes (GTs 1 ~ 74) (Fig. 2), 57 of which were represented by unique strains, and the remaining 17 were MLVA-16 genotypes, being shared between two and four isolates (Fig. 2). The 17 shared genotypes comprised 44 clustered strains for a clustering rate of 43.00% (44/101; Fig. 2). The results indicated a co-existing transmission pattern featuring sporadic cases and limited epidemiologically related cases. Ten shared genotypes (GT1, 2, 3, 13, 30, 47, 52, 68, 69, and 70) (Fig. 2, black shadow) contained strains from the same location and isolation time, implying that multiple outbreak events were caused by a common source (Fig. 2). Other shared GTs consisted of isolates from different regions. Further investigations are needed to reveal the pattern of transmission of Brucella strains in this province.

Fig. 2

The dendrogram of the MLVA-16 of 101 Brucella melitensis strains in Sichuan Province. Key: strains code; GT: MLVA-16 genotypes; MLVA-8: MLVA-8 genotypes; MLVA-11: MLVA-11 genotypes; Year: the isolation dates of strains; isolated in: the locations of isolated strains

Genetic comparison of Brucella strains from China and the global scale

An MLVA-16 comparison analysis was conducted to investigate the molecular epidemiological relationships among strains at the national scale. The MLVA comparison of 801 B. melitensis identified 36 shared MLVA-16 genotypes (1–36) (Table S4) between strains from Sichuan and 22 other provinces (Fig. 3). A total of 51 strains (50.5%, 51/101) from the present study shared MLVA-16 genotypes with strains from other provinces (cities). The provinces that shared strains with MLVA-16 genotypes included Ningxia (n = 40), Qinghai (n = 36), and Inner Mongolia (n = 20) (Fig. 3), followed by Shandong (n = 7), Shanxi (n = 6), Shaanxi (n = 5), and Xinjiang (n = 5). These provinces were historic endemic areas of human brucellosis, and parts of each are contiguous with Sichuan Province. On a global scale, multiple MLVA-16 genotypes were shared by strains from the present study, and strains from India (n = 1), Mongolia (n = 2), and Kazakhstan (n = 14) (Fig. 4). Strains from the present study exhibited high genetic similarity with strains from many countries in Asia, including Saudi Arabia, Turkey, Kuwait, Afghanistan, Lebanon, Iraq, Syria, and Israel (Fig. S5).

Fig. 3

A minimum spanning tree of the genetic relationships of B. melitensis strains from a national survey constructed using the PHYLOVIZ 2.0 software. Note: The 51 B. melitensis strains in this study formed 36 identical MLVA-16 genotypes with strains from 22 different provinces in China. The numbers in circles refer to the MLVA-16 genotypes, and the numbers in lines indicate the differences in loci between two adjacent MLVA-16 genotypes. Strains from different provinces (cities) are marked with colors. The strains from the present study are marked with red, with light pink for Ningxia, blue for Qinghai, light green for Inner Mongolia, purple for Shandong, blue for Shanxi, dark green for Shaanxi, and deep blue for Xinjiang

Fig. 4

Minimum spanning tree of genetic correlations of B. melitensis strains at a global scale. Note: The MST was constructed based on the MLVA-16 data from 1296 B. melitensis strains using the Bionumerics software 8.0. The strains from the present study are marked with yellow and highlighted with yellow circles. The strains from other countries coded with different colors such as light yellow for Mongolia, light green for Kazakhstan, and indigo for India

Discussion

In this study, bacteriology, bio-typing, and two genotyping approaches were applied to discriminate 101 Brucella strains isolated from Sichuan Province. At present, B. meltenesis bv. 3 is the dominant circulating species in Sichuan Province, and this strain has spread to 16 out of 21 cities. Previous studies demonstrated that B. melitenis has continuously spread from north to south, posing a threat to public health [6]. Ruminants are the preferred reservoirs, implying that imported infected sheep and goats are the main source of human brucellosis in the region. Therefore, control measures for brucellosis are urgently needed to contain the spread of the pathogen. This conclusion is consistent with the epidemic situation in Sichuan Province, where human brucellosis is a re-emerging disease. Sporadic cases have occurred since 2014, with reported cases increasing from 25 in 2014 to 417 in 2023. In addition, Sichuan Province borders many provinces from the west, including Shaanxi, Qinghai, Gansu, and Yunnan provinces, where human brucellosis has a high prevalence [18]. Also, the breeding numbers of sheep and goats are high in Sichuan owing to the development of animal husbandry to increase the incomes of local farmers and herdsmen.

The sheep stock increased from 1689.19 (tens of thousands) in 2013 to 1761.13 (tens of thousands) in 2016, and then declined to 1529.9 (tens of thousands) in 2022. Mutton and lamb production increased from 22.45 (ten thousand tons) in 2013 to 27.10 (ten thousand tons) in 2023 (data from the Bureau of Statistics, China). Due to the development of the breeding industry, the local population experienced indirect and direct contact with infected ruminants, and this has driven the re-emergence of brucellosis and local epidemics.

Strains belonging to ST8 are distributed in northern and northeast China, including Inner Mongolia, Xinjiang, Gansu, Qinghai, Shaanxi, and Heilongjiang provinces [19]. The latest reports show that B. melitensis ST39 was first isolated in Zunyi City, Guizhou Province in 2013, and many strains belonging to the ST group have been isolated in this province [20]. The isolated B. melitensis ST8 and ST39 in Sichuan imply that the source of infection was from outside the province; however, the details concerning the epidemiological links among these strains need further discrimination by WGS-SNP [21]. Moreover, two new STs were detected in this study, and each represented only one strain. Sequence comparison between ST8 and ST101 showed that they differed at the glk locus, while ST118 differed at the trpE locus. The source and epidemic status of the two new sequence types require further investigation.

The predominance of the 116 and 125 of MLVA-11 genotypes of Brucella in Sichuan Province was consistent with the pattern in most provinces of China. Genotype 116 is the predominant genotype of B. melitensis that continues to expand from northern to southern China [6]. MLVA-16 cluster analysis revealed a co-existing transmission pattern dominated by sporadic cases and multiple outbreak events caused by a common source. This result was consistent with the epidemic status of human brucellosis in this province, where the incidence rate was 0.498/100,000 in 2023. In contrast, MLVA genotyping confirmed that the epidemic pattern was dominated by a multipoint outbreak in Ulanqab, Inner Mongolia [22]. This result implied that human brucellosis is in the early epidemic stages in this region, and thus timely and targeted control measures are needed. However, the nationwide MLVA comparison analysis revealed that 51 strains (50.5%, 51/101) in this study involved 36 shared MLVA-16 genotypes with strains from 22 other provinces. These data imply that the dominant MLVA genotypes of B. melitensis have circulated among various provinces and that the shared genotype of MLVA-16 of B. melitensis in Sichuan province and other places in China confirmed the inference of extra-provincial importation, indicating regional transmission of B. melitensis strains from the same lineages that is possibly associated with the livestock trade [23]. This indicates that massive-scale farming may be increasing the risk of transmission of brucellosis, and that most local farmers may not be aware that introducing infected ruminants will spread the disease. Furthermore, matching of local genotypes with global genotypes revealed that some local isolates had an East Mediterranean origin, and strains in this study displayed a level of high genetic homogeneity with strains from other countries in Asia. The East Mediterranean lineage of B. melitensis is predominant in Asia [24]. B. melitensis populations represent a public health challenge in countries along with Silk Road [23]. And also, B. melitensis biovar 3 is the most prevalent biovar in Middle Eastern and Asian countries, which caused the most human cases [25].

Many studied underscored that necessity of whole genome sequencing (WGS) to gather more information for improved epidemiological investigations [26, 27]. A core-genome single-nucleotide polymorphism (cgSNP) analysis revealed a potential connection between the Turkish isolates of B. melitensis and those from Sweden, Israel, Syria, Austria, and India [28]. Mongolian B. melitensis isolates had high genetic similarity to Chinese strains, likely due to the geographical proximity [29]. Extensive genotype similarities were observed between strains from Inner Mongolia, Kazakhstan, Mongolia, and Turkey. Brucellosis is a significant zoonotic infection in ruminants in Pakistan. In a survey, the results point toward a persistent zoonotic threat in cattle in the district and the potential spillover infection when close contact with other animal species [30]. These countries were key members of the grassland silk road, and long-term trade in small ruminants (sheep) in these countries has possibly promoted the spread of Brucella spp. in these regions [31]. This explains the genetic links between Brucella strains from neighboring regions/countries and the isolates from Sichuan Province, China. Therefore, animal quarantines should be strengthened to prevent the spread of Brucella to adjacent countries and regions [32]. Brucellosis in the Nile River Basin countries (Egypt, Sudan, Ethiopia, and Tanzania) is highly prevalent and endemic. There are several factors behind the failure of eradication of Brucella in these countries, including the lack of cooperation between policymakers, health officials, veterinary sectors, and farmers; and inappropriate disposal of aborted materials, insufficient government compensation for infected animals, and public vaccination reluctance, open borders and uncontrolled animal movements, lack of proper diagnostics and bio-typing methods and tools [33]. Therefore, implementing a strict testing and slaughter policy is needed to eliminate the infected ruminants and reduce the incidence of human brucellosis [34]. In addition, the movement of infected animals should be controlled, where financial and material investment by government departments should aid farmers’ participation in control programs.

Conclusion

Although the Brucella isolates in Sichuan Province were all B. melitensis bv. 3, the sequence types and MLVA-16 genotypes indicated high genetic diversity, and the present epidemic of human brucellosis is dominated by sporadic cases. The control of the disease is challenging due to the extensive MLVA-16 genotype comprising strains from Sichuan and 22 other provinces. Most of the shared strains were from the northern region, where brucellosis occurs at a high prevalence. A WGS-SNP-based analysis employed to investigate strains at a regional and national scale will provide data for a better understanding of the re-emergence of human brucellosis in Sichuan and other regions. The control of human brucellosis in Sichuan Province requires strict, multi-sectoral cooperation to limit outbreaks of brucellosis in domestic animals.