Enzootic situation and molecular epidemiology of Brucella in livestock from 2011 to 2015 in Qingyang, China

A large-scale survey was conducted in domestic animal populations from 2011 to 2015 in Qingyang, China. A total of 448,398 animals from different districts of Qingyang were tested for the presence of Brucella-specific antibodies using the Rose Bengal Plate Test (RBPT) and the Standard Agglutination Test (SAT). From 2011 to 2015, the yearly average positive rates were between 0.04 and 4.75% in the eight counties tested. In addition, the prevalence rates were between 0 and 9.96% in these eight counties. Sheep was the dominant host of Brucella in Qingyang, and the prevalence rate in sheep (2.74%) was higher than those in the other animals tested. Identification of 10 Brucella isolates from sheep confirmed that the epidemic strains were B. melitensis biovar 3 (n = 9) and B. melitensis biovar 1 (n = 1). MLVA-11 (multiple-locus variable-number tandem repeat analysis) analysis of the 10 isolates showed three genotypes: genotype 116 (n = 8), genotype 115 (n = 1) and genotype 136 (n = 1). Furthermore, analysis of the whole-genome sequences of the representative B. melitensis strain QY1 indicated that this isolate was closely related to isolates from China and India. The results of serum epidemiology confirmed that the region of northern Qingyang was a critical Brucella epidemic area and that the disease showed a rising trend, especially from 2013 to 2015. An analysis of the isolate genotypes suggested that sheep brucellosis mainly resulted from conventional B. melitensis (East Mediterranean group), although the external strain (American group) also occurred in Qingyang.


Introduction
Brucellosis, caused by Brucella spp., is a zoonotic bacterial disease that has serious implications for human and animal health 1 . It is a highly contagious disease that affects livestock such as cattle, sheep and pigs. The symptoms include abortion, infertility, decreased production and lameness in animals and undulating fever with arthralgia in humans 2 . Recently, this disease has reemerged in some regions of China [3][4][5][6] . From an epidemiological perspective, the major cases in animals have been sheep infected by B. melitensis. Human brucellosis has re-occurred in the past several years 7,8 . In Qingyang, four cases were detected from 2006 to 2008, while from 2009 to 2012, the number of positively detected Brucella cases was between 25 and 57. In 2013 and 2014, however, the cases of human brucellosis numbered 161 and 457, respectively 9 . Among animal cases, from 2009 to 2013, the seroprevalence was <1%, although brucellosis was already occurring in domestic animals 10 .
Over the past decade, the number of breeding stock of both cattle and swine were approximately 400,000 in Qingyang. However, sheep farming was developed vigorously with main farmers feed. From 2007 to 2013, more than two million sheep were bred, and the number of sheep breeding stock was over 4 million in both 2014 and 2015 in Qingyang. From the beginning of this century to 2016, although animal brucellosis occurred in an increasing number of cases, the local control strategy was to cull the infected animals that presented as serologically Brucella positive. After 2016, vaccination was required in the area because of changes in government control policy, which were based on an assessment of the development of brucellosis in animals. Therefore, detailed knowledge of the epidemiological situation has become very important for the vaccination program and for assessment of effective prevention. This study was aimed at understanding the enzootic situation of Brucella in livestock and the characteristics of epidemics. Comprehensive Brucella seroepidemiological surveys were performed in livestock epidemics of brucellosis from 2011 to 2015 in Qingyang,  China. To further illuminate the molecular features of epidemic Brucella, bacterial isolation and whole-genome sequencing were also performed.

Domestic brucellosis in Qingyang from 2011 to 2015
Of a total of 448,398 serum samples from different animals, 11,654 (2.6%) were anti-Brucella positive. The yearly positive rates in swine, yellow cattle (native breed), dairy cows (Holstein), sheep and dogs were 0%, 0.04%, 1.66%, 2.74% and 0.14%, respectively. Sheep and dairy cows with Brucella had higher infection trends, although the positive rates in dairy cows dipped to a downward trend between 2013 and 2014. Compared with that in sheep and dairy cows, brucellosis in yellow cattle and dogs maintained sporadic onset, and no cases were found in swine (Table 1). From 2011 to 2015, the average positive rates in each year were from 0.04 to 4.75% in eight counties. In addition, the yearly prevalence rates were among 0 to 9.96% in these eight counties, and the disease emerged with an obvious rising trend, especially in the years 2013 and 2015. However, brucellosis in animals showed regional differences (Figs. 1a, b). Overall, domestic brucellosis in Qingyang predominantly infected sheep. Huanxian, located in the north Qingyang district, was a serious animal Brucella epidemic area. In the years 2013 and 2014, 6 (0.04%) of the 13,793 total yellow cattle were brucellosis positive, and 0.14% (1/698) seroprevalence was detected in dogs (Table 1).

Dairy cow brucellosis in Qingyang
Serological examination revealed that dairy cow brucellosis mainly emerged in the southern region of Qingyang, including Ningxian, Zhenyuan and Zhengning counties, whereas no dairy cow Brucella cases were found in the northern region of Qingyang. The yearly prevalence rates of Brucella-infected dairy cows was from 0.89 to 4.86% from 2011 to 2015 without incidence trends, and the highest positive detection rate was 20.12% in Zhenyuan in 2015 (Figs. 1c, d).

Sheep brucellosis in Qingyang
To clarify the situation of Brucella in sheep in Qingyang district, an analysis of the epidemic in eight counties showed that the Brucella prevalence in each year ranged from 0 to 10.21% (Fig. 2), and the average was found to be 0.33 to 4.79% from 2011 to 2015 (Figs. 2a, b). Furthermore, the prevalence of Brucella showed an obvious rising trend, and the disease displayed a significant regional difference (Figs. 2c, d). According to the geographical locations of the eight counties, Brucella in sheep showed an increasing gradient from south to north. Three southern counties, Zhengning, Ningxian and Xifeng, had lower positive rates (<2%) in sheep, but the positive rate in northern Huanxian was the highest (9.7%). In the other counties, including Zhengyuan, Qingcheng, Huachi and Heshui, the positive rates of brucellosis were from 2.65 to 4.17% (Fig. 3). Although the average seroprevalence in Huanxian was significantly higher than those in other counties, no regional difference (P > 0.05) was displayed because of the outbreak in 2015.

Isolation and identification of Brucella strains in sheep
From 55 spleens of aborted samples, a total of 10 isolates were obtained. Specifically, five and two isolates were isolated from Huanxian and Zhenyuan, and one isolate was isolated from each of Huachi, Heshui and Zhengning. Biochemical testing and AMOS-PCR showed that nine biovars were identified, of which one was B. melitensis biovar 1 and nine were B. melitensis biovar 3 (n = 1; Fig. 4).

Phylogenetic and evolutionary relationships of B. melitensis QY1
Analysis of the whole-genome sequences showed that B. melitensis QY1 has two circular chromosomes with 2,125,648 bp in chromosome I and 1,185,604 bp in chromosome II. The results of a phylogenetic tree based on the whole genome demonstrated that B. melitensis QY1 is close to B. melitensis M28, B. melitensis NI and B. melitensis M5-90, which originated in China, as well as B. melitensis bv. 3 str. Ether, a strain that was isolated from India (Fig. 6).

Discussion
Qingyang is located at the intersection of Gansu, Shaanxi and Ningxia provinces and is near Inner Mongolia. The latter two provinces have very large livestock populations with high prevalences of Brucella 3,11 . In the present study, we used serological methods to investigate animal brucellosis from 2011 to 2015 in Qingyang, China. The results of seroprevalence assays confirmed that Brucella spp. was mainly epidemic in sheep and dairy cows, whereas other animals had a low incidence. The results also showed that the disease has attained an obvious rising trend. Similarly, the number of human cases of brucellosis has also increased 9,12,13 . These results imply that the prevalence of brucellosis in animals, especially in sheep, has caused an increase in the human disease. In contrast, although Brucella infected a large number of sheep and dairy cows, yellow cattle presented very few cases.
From 2011 to 2015, brucellosis in dairy cows emerged without incidence trends in Qingyang. There are two possible explanations for this phenomenon. The breeding stock of dairy cows was <3000, so the disease was easy to control. In addition, the local control strategy was to cull infected animals, which largely prevented the expansion of Brucella in the short term. From an epidemiological perspective, sheep brucellosis emerged broadly in Qingyang, and dairy cow brucellosis occurred as outbreaks in some dairy farms. It was implied that brucellosis in dairy cows came from infected sheep.
Currently, the number of sheep breeding stock is over 4 million in Qingyang, and almost half of the sheep were raised in Huanxian. From 2011 to 2014, the seroprevalence of brucellosis in Huanxian was <2%, but an outbreak of the disease emerged in the year 2015. The most important factors may be convenient transportation and the prevalence of Brucella in most regions of China, as well as monitoring of animals with Brucella. Meanwhile, in the past few years, the prevalence of brucellosis in animals and humans has increased each year in Ningxia province, which is adjacent to northern Huanxian 11 . This result may be explained by the fact that frequent livestock exchange occurs between these areas. The outbreak of brucellosis in 2015 further reflected a lack of control in the movement of infected animals between the regions and the supervision of the market. This hypothesis suggests that inspection measures to ensure the quarantine of livestock exported from endemic regions should be strengthened 7 .
To clarify the molecular features of sheep Brucella in Qingyang, 10 isolates were revived from Huanxian, Zhenyuan, Huachi, Heshui and Zhengning. Genotyping by MLVA-11 showed that the isolates were divided into three relatively independent clusters (cluster A, B and C; Fig. 5). Cluster A, containing five isolates, was very close to major strains that were responsible for the dominant form of sheep brucellosis in China. However, the isolates in Cluster B were closely related to a few strains isolated from Ningxia, Inner Mongolia, Heilongjiang and Shandong 8 . An isolate in Cluster C belongs to B. melitensis biovar 1, which was much closer to Brucella strains isolated from Hainan, China, in 2013 8 . The results of MLVA-11 confirmed that major isolates were genotype 116 (n = 8), a common genotype derived from Brucella strains in the East Mediterranean. Moreover, this genotype was also responsible for the vast majority of Brucella infections in China. However, an isolate (genotype 136) from an American lineage was identified and observed in this study. Based on the analysis of these results, it is suggested that sheep brucellosis in Qingyang occurred not only due to a primary epidemic focus but also because of invasion by external bacterial strains with the introduction of sheep.
B. melitensis QY1, as a representative strain, was sequenced and analyzed. Sheep brucellosis re-emerged in Qingyang with rising trends. Therefore, a wide range of vaccination strategies were adopted starting in 2016. China has applied a vaccination program using live attenuated strains S2 (B. suis biovar 1 S2) and M5 (B. melitensis M5-90) [14][15][16] , but the problem is to distinguish positive sera due to vaccination from positive sera due to field strain infection. Modified live attenuated vaccines have been developed using differential diagnosis, but so far there are no commercial applications. Studies of molecular epidemiology and strain characteristics can provide support for vaccine selection and development.

Ethics statement
This study involves the investigation and isolation of Brucella spp. using serology and modern typing methods. All animal experiments in the study were approved by the Animal Ethics Committee of Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences, and were conducted strictly according to the guidelines.  (Table 1). These samples were collected from eight counties, including Xifeng (32,478), Zhengning (30,904), Ningxian (56,228), Heshui (62,906), Zhenyuan (44,576), Qingcheng (48,212), Fig. 6 Phylogenetic tree based on the whole genomes of nine selected B. melitensis strains. The evolutionary history was inferred using the neighbor-joining method. The optimal tree (with branch length sum = 0.00316668) is shown. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree Huachi (75,396) and Huanxian (98,213) (Fig. 1, Table 1 and Table S2). The sample collection rates in sheep were from 2 to 5% (proportion of total population), the dairy cows were fully sampled, and sampling rates in the other livestock were higher than 20%. All sample collection was performed to cover the full area of Qingyang district, which was sampled randomly. All serum samples were transported to the laboratory, where they were stored at −20°C until processing and detection by Rose Bengal Plate Test (RBPT), and the positive serum samples were reconfirmed by a Standard Agglutination Test (SAT) 17 (all the reagents were purchased from China Institute of Veterinary Drug Control).

Isolation of Brucella
To study the molecular epidemiological characters, 55 spleens of aborted sheep were collected at random from eight counties and subjected to a brucellosis purification progress. Spleens were crushed and cultured on Brucella serum dextrose agar composed of Brucella medium base (supplemented with Brucella selective antibiotic, OXOID, England), and 5-10% heat-inactivated horse serum. Plates were incubated with and without 5-10% carbon dioxide at 37°C after inoculation with sample materials. The plates were examined after 3-15 days for bacterial growth. A single clone was chosen for identification.

Identification of isolates
The obtained single bacterial clones were identified by biochemical testing according to the standard strain identification method [17][18][19] . The carbon dioxide (CO 2 ) requirement was tested on Brucella serum dextrose agar with and without CO 2 during the first isolation. Agglutination by A, M and R monospecific antisera were detected by mixing the antisera with the isolate after dilution of the colony (all reagents were purchased from The Chinese Center for Disease Control and Prevention). Total genomic DNA was extracted using a DNeasy Blood & Tissue Kit (Qiagen, Germany) according to the manufacturer's instructions. The DNA extracted from all isolates was stored at −20°C. AMOS-PCR was performed as previously described 20,21 .

Molecular epidemiological tracking of isolates by MLVA-16
MLVA was performed as described previously [22][23][24] . Sixteen primer pairs were included in panel 1 (bruce06, bruce08, bruce11, bruce12, bruce42, bruce43, bruce45 and bruce55), panel 2A (bruce18, bruce19 and bruce21), and panel 2B (bruce04, bruce07, bruce09, bruce16 and bruce30). Polymerase chain reaction (PCR) amplification was performed in a 50 μL reaction system, containing 1 μL DNA template, 1 μL of each primer, 25 μL of 2 × Premix Taq, and 22 μL of ddH 2 O. The reaction conditions were as follows: initial denaturation at 94°C for 5 min; followed by 40 cycles of 94°C for 30 s, 60°C for 30 s and 72°C for 1 min. At last, the samples were incubated for an additional 10 min at 72°C and were then stored at 4°C. Five microliters of each PCR product from panel 1 was loaded onto a 2% agarose gel with 0.5% ethidium bromide, whereas those from panels 2A and 2B were loaded into 3% agarose gels also containing 0.5% ethidium bromide. The bands were visualized under UV light and photographed. The PCR products for the 16 loci were denatured and resolved by capillary electrophoresis on an ABI Prism 3130 automated fluorescent capillary DNA sequencer (Applied Biosystems). Fragments were sized following comparison with a ROX (carboxy-X-rhodamine)-labeled molecular ladder (MapMarker 1000; BioVentures Inc., Murfreesboro, TN, USA) and Gene Mapper software version 4.0 (Applied Biosystems). Band intensities were estimated using BioNumerics version 7.6 (Applied Maths, Belgium), and repeat units were then obtained according to the published allele numbering system 22,23 . Clustering analysis was performed using the same software based on the categorical coefficient and unweighted pair group using arithmetic averages method. The resulting genotypes were compared using the web-based Brucella 2016 MLVA database (http://mlva.u-psud.fr/) and published papers 23 . Advanced clustering analysis was based on the MLVA for categorical data in the BioNumerics system.

Genomic sequencing and evolutionary relationships of B. melitensis QY1
A dominant epidemic isolate, B. melitensis QY1, was used for whole-genome sequencing. Genomic DNA was extracted with the sodium-dodecyl sulphate (SDS) method. The harvested DNA was detected by agarose gel electrophoresis and quantified by Qubit (Life Qubit ® Fluorescence Quantitative Instrument, USA). The genome of B. melitensis QY1 was sequenced by singlemolecule real-time (SMRT) technology. Sequencing was performed at Beijing Novogene Bioinformatics Technology Co., Ltd. The low-quality reads were filtered by SMRT 2.3.0 25,26 , and the filtered reads were assembled to generate one contig without gaps (accession nos. CP022204.1-CP022205.1) 27 . The evolutionary relationships were evaluated in comparison with B. melitensis M28, B. melitensis M5-90, B. melitensis NI, B. melitensis bv. 1 str. 16M, B. melitensis bv. 3 str. Ether, B. melitensis biovar abortus 2308 and B. melitensis ATCC 23457 28 (Table S3). Based on the whole genomes of nine selected B. melitensis strains, the amino-acid sequences of core genes were analyzed by the CD-HIT software for rapid clustering of similar proteins with a threshold of 50% pairwise identity and 0.7 length difference cutoff. A phylogenetic tree was constructed from the orthologous gene sequences by TreeBeST using the neighbor-joining method, and the number of bootstrap samples was 1000 29,30 .

Statistical analysis
The data were transferred to Microsoft Excel and evaluated using a two-tailed unpaired t-test to compare the two groups; P-values <0.05 were considered to be significantly different.