Genomic characters of Anaplasma bovis and genetic diversity in China

ABSTRACT

The emergence of Anaplasma bovis or A. bovis-like infection in humans from China and the United States of America has raised concern about the public health importance of this pathogen. Although A. bovis has been detected in a wide range of ticks and mammals in the world, no genome of the pathogen is available up to now, which has prohibited us from better understanding the genetic basis for its pathogenicity. Here we describe an A. bovis genome from metagenomic sequencing of an infected goat in China. Anaplasma bovis had the smallest genome of the genus Anaplasma, and relatively lower GC content. Phylogenetic analysis of single-copy orthologue sequence showed that A. bovis was closely related to A. platys and A. phagocytophilum, but relatively far from intraerythrocytic Anaplasma species. Anaplasma bovis had 116 unique orthogroups and lacked 51 orthogroups in comparison to other Anaplasma species. The virulence factors of A. bovis were significantly less than those of A. phagocytophilum, suggesting less pathogenicity of A. bovis. When tested by specific PCR assays, A. bovis was detected in 23 of 29 goats, with an infection rate up to 79.3% (95% CI: 64.6% ∼94.1%). The phylogenetic analyses based on partial 16S rRNA, gltA and groEL genes indicated that A. bovis had high genetic diversity. The findings of this study lay a foundation for further understanding of the biological characteristics and genetic diversity of A. bovis, and will facilitate the formulation of prevention and control strategies.

KEYWORDS:

• Anaplasma bovis
• tick-borne disease
• whole-genome sequencing
• genetic diversity
• goat
• China

Introduction

Anaplasma species in the genus Anaplasma of the family Anaplasmataceae, order Rickettsiales are Gram-negative obligate intracellular Alphaproteobacteria mainly transmitted by ticks [1,2]. Six species are previously known in the genus Anaplasma, including A. bovis, A. centrale, A. marginale, A. ovis, A. phagocytophilum and A. platys. Two newly recognized species, A. capra and A. odocoilei [3,4], and a number of new or unclassified candidate species have emerged in the past two decades [5]. The Anaplasma species infect a wide range of animal hosts, and more and more human infections with various Anaplasma species have been continuously reported.

Anaplasma bovis, previously called Ehrlichia bovis, was initially detected in cattle in Brazil in 1936. It is widely distributed in Asia, Africa, America and southern Europe, and can cause cattle disease with the manifestations of fever, weight loss, lower milk production, enlargement of lymph nodes, and occasionally abortion or death [1]. Anaplasma bovis has been detected in domestic and wild animals, including goats, sheep, deer, dogs, cats, rabbits, Mongolian gazelles, raccoons, monkeys, and small mammals [6–12]. The identification of A. bovis in cynomolgus monkeys from Malaysia suggests its potential infectivity to humans [12]. The recent reports of human infections with A. bovis in China [13,14] and human A. bovis-like infections in the United States of America [15] have raised concerns about its public health importance.

With the development of sequencing technology, a total of 24 Amarginale genomes [16,17], 32 A. phagocytophilum genomes [18,19], two A. ovis genomes [20], two A. capra [21], an A. centrale genome [22], and an A. platys genome [23] have been sequenced and deposited in GenBank (https://www.ncbi.nlm.nih.gov/genome). In spite of its wide-distribution, extensive detection in ticks and animals around the world, and confirmed human infections, there is no A. bovis genome up to now. In this study, we assembled a whole genome of A. bovis from a blood sample of an infected goat in Shandong Province, China. We subsequently characterized the genomic features of A. bovis in comparison to other species in the genus Anaplasma and investigated the genetic diversity of A. bovis in China.

Materials and methods

Sample collection and DNA extraction

EDTA-anticoagulated blood samples were collected from free-range goats at Yantai City of Shandong Province, China in 2022. Total DNA was extracted from each goat blood sample using a blood DNA extraction kit (Omega Biotek Inc., Norcross, GA, USA) according to the manufacturer’s instructions.

PCR assays and product sequencing

A nested PCR targeting the Anaplasma 16S rRNA gene [24] was performed for the blood samples to screen Anaplasma infections among the goats. All the positive blood samples were subsequently tested using a specific PCR assay targeting the citrate synthase (gltA) and chaperonin GroEL (groEL) genes [7,25] to investigate the genetic diversity of A. bovis (Table S1). All the PCR products were sequenced by Sanger sequencing to confirm the amplification results.

Metagenomic sequencing of a goat blood sample

A goat blood sample positive for 16S rRNA, gltA and groEL genes of A. bovis was used for metagenomic sequencing. Total DNA was extracted from the positive goat blood and purified to satisfy the quantity and quality of sequencing. The sequencing library was then constructed following the whole genome sequencing library preparation protocol. The paired-end (PE) library was sequenced with a read length of 2 × 150 bp on a DNBseq-T7 platform at Grandomics Gene Technology Beijing Co. Ltd (Beijing, China).

Assembly and annotation of A. bovis genome

Cleaned sequence reads were mapped to the genome of goat (Capra hircus) (GenBank accession number GCF_001704415) by bowtie2 v2.4.1 [26] with default options for eliminating the host genome. The remaining reads which were not matched to the goat genome were converted to BAM files by SAMtools v1.9 [27]. Host-free reads were inputted to SPAdes v3.15.2 [28] with -a meta parameter to assemble reads to contigs and scaffolds. MetaBAT2 v2:2.15 [29] was used for contig binning. Prokka v1.13.3 [30] and check M v1.1.3 [31] with lineage_wf argument were used for genome annotation and genome completeness assessment. Orthofinder v2.5.4 [32] was utilized to find single-copy orthologue sequences and orthogroups. The local pair mode of software mafft v7.505 was used for multi-sequence alignment [33], Gblocks v0.91b for divergent region elimination [34], and iqtree v2.2.0.3 for building phylogenetic tree [35]. UpsetR was used to visualize share, unique and intersection situations in orthogroups [36]. KOfamKOALA v1.4.0 and eggNOG-mapper v2.1.7 were used for function annotation of selected orthogroups and across the whole genome [37,38].

Analysis of genetic diversity of A. bovis

Information on the collection location and host of each A. bovis sequence in China was extracted from GenBank and its related publication. Longitude and latitude were transformed from the location by API of Gaode map (https://lbs.amap.com/tools/picker). Complemented sequences were aligned by mafft v7.505 local pair mode [33], and phylogenetic trees were made by iqtree v2.2.0.3 [35]. The record in China was visualized, and the phylogenetic tree was displayed using the R package sf and ggtree [39,40], respectively.

Results

Assembly of A. bovis from an infected goat blood

A goat blood sample positive for 16S rRNA, gltA, groEL genes of A. bovis by PCR amplification (Table S1) was chosen for next-generation sequencing. Total DNA was extracted from 5 ml blood, and subsequently processed according to the whole-genome sequencing library preparation protocol (DNBSEQ) to build a sequencing library. The NGS generated 78 million 150-bp Illumina reads from the sample. About 97% of total reads were mapped to the goat genome and discarded. The remaining reads were then de novo assembled into contigs and scaffolds using the SPAdes 3.15.3 with meta parameters [27,28]. After assembly and binning, a genome of A. bovis was identified by CheckM, and named A. bovis str. BIME (GenBank accession no. JAWLNY000000000), because it was initially obtained by the research team at the Beijing Institute of Microbiology and Epidemiology, China.

Genomic characteristics and phylogenetic position of A. bovis str. BIME

The genome size of A. bovis str. BIME was 1.05 Mbp with 95.30% completeness, which was the smallest among those in the genus Anaplasma, including A. capra, A. centrale, A. marginale, A. ovis, A. phagocytophilum and A. platys. GC content (42.94%) of this genome was the second lowest in the genus Anaplasma, which was only higher than A. phagocytophilum. There was no correlation between genome size and GC content of different species (Figure 1(A)). According to prokka annotation, there were 923 genes in the genome of A. bovis str. BIME, 844 of which were coding sequences. The genome had 36 tRNAs, one tmRNA and 42 pseudogenes (Table 1). It showed that A. bovis possessed the smallest genome, second smallest GC content and pseudogene counts in the genus Anaplasma. A positive correlation was observed between genome size and number of pseudogenes (Figure S1).

Figure 1. Basic information of A. bovis genome. (A) Scatter plot of Anaplasma representative genome. (B) Circos plot of A. bovis genome. From inner circle to outer circle representative GC content, GC skew, sequence depth, proteins of + strand, proteins of – strand. Red arrows symbolize virulence factor. (C) Phylogenetic tree of ML method and 1000 replications based on single copy orthologue of all Anaplasma genome.

Table 1. Basic information for representative strain of Anaplasma species.

Scientific name	A. bovis str. BIME	A. capra str. BIME2	A. centrale str. Israel	A. ovis str. Haibei	A. marginale str. Florida	A. platys str. S3	A. phagocytophilum str. JM
Accession number	JAWLNY000000000	GCF_025628805.1	GCF_000024505.1	GCF_002214625.1	GCF_000020305.1	GCF_012790675.1	GCF_000439775.1
Genome size (bp)	1,052,900	1,059,758	1,206,806	1,214,674	1,202,435	1,196,811	1,481,598
GC content (%)	42.94	48.29	49.98	48.89	49.77	45.18	41.64
L50	19	5	1	1	1	1	1
N50 (bp)	16,948	100,463	1,206,806	1,214,674	1,202,435	1,196,811	1,481,598
Genes	923	930	1059	1066	1056	1054	1276
CDS	844	859	906	868	864	929	823
tRNA	36	37	37	37	37	34	37
tmRNA	1	–	1	1	1	1	1
rRNA	–	3	3	3	3	3	3
Pseudogene	42	31	112	157	151	87	412
Completeness	95.30%	–	–	–	–	–	–

The Circular map depicted the basic characteristics of A. bovis genome (Figure 1(B)). From the inner to the outer side, circles displayed GC content, GC skew, sequencing depth, coding sequences on the reverse and forward strands, and genome length. These elements were evenly distributed in this genome with no huge gaps and saltation. A phylogenetic tree based on single copy orthologue sequences showed that A. bovis was distinctive from other Anaplasma species, and closer with A. platys and A. phagocytophilum than intraerythrocytic Anaplasma species, including A. capra, A. centrale, A. marginale and A. ovis (Figure 1(C)).

Annotation of A. bovis genome str. BIME in comparison to other Anaplasma species

To clarify the difference among Anaplasma species, eggNOG-mapper was used to annotate proteins. A bar plot, which counted proteins of representative genome classified by COG category showed that the top three categories of A. bovis were J (Translation, ribosomal structure, and biogenesis), S (Function unknown) and C (Energy production and conversion), same as other Anaplasma species. Notably, the number of proteins in A. bovis was fewer than those in other Anaplasma species among most categories. Anaplasma bovis and A. platys had much fewer proteins in category M (Cell wall/membrane/envelope biogenesis) compared to other species, suggesting lower functions of the two species in cell membrane biogenesis (Figure 2(A), Table S2).

Figure 2. Functional annotation of A. bovis genome. (A) COG annotation of Anaplasma genome. (B) UpsetR plot of intersections between Anaplasma representative genome. Circles which were connected means these species of Anaplasma sharing orthogroups. (C) Heatmap of Anaplasma representative genomes virulence factor holding situation. Size of circle representative coverage of virulence factors, colour of circle representative identity of virulence factors.

To further explore differences among representative Anaplasma genomes, Orthofinder was used to find orthogroups, and UpsetR was used for visualization. All the Anaplasma genomes shared 651 orthogroups. Notably, A. bovis missed 51 homologous groups compared to the other six Anaplasma species. Anaplasma bovis, A. platys, and A. phagocytophilum lacked 77 homologous groups in comparison to the erythrocytic Anaplasma species, including A. capra, A. centrale, A. ovis, and A. marginale, but had three specific orthogroups. Anaplasma bovis possessed 116 unique orthogroups (Figure 2(B)). The functions of these 116 unique orthogroups were further annotated by using KEGG, and they mainly focused on genetic information processing, Protein families: metabolism and Protein families: signalling and cellular processes (Figure S2).

The bubble plot shows the differences in amino acid sequences and coverage among Asp14, OmpA, VirB2-7, VirB2-8, VirB3, VirB4, VirB6, VirB7, VirB8, VirB9, VirB10, VirB11, VirD4 and AnkA between A. bovis and the erythrocytic Anaplasma species (A. centrale, A. marginale, A. ovis, and A. capra). Anaplasma bovis was most similar to A. platys in most virulence factors except for VirB6-3 and VirB6-4. Compared with A. phagocytophilum, the coverage and similarity of all virulence factors in A. bovis were significantly lower, suggesting that A. bovis might be less pathogenic to humans and animals (Figure 2(C)).

Diversity and prevalence of A. bovis in China

In this study, the specific primers of 16S rRNA, gltA and groEL were used for PCR detection, and Sanger sequencing was performed to screen for A. bovis from 29 goats in Shandong Province. A total of 23 were positive for 16S rRNA, 27 for gltA and 6 for groEL (Table S3). To ensure the reliability of test results, a sample was considered positive only if two genes in the sample were amplified to be positive. As a result, the positive rate was 79.31% (95% CI: 64.6% ∼94.1%). The sequencing results of 16S rRNA, gltA and groEL were deposited to GenBank with accession numbers OR717419-OR717441 for 16S rRNA, OR731334-OR731360 for gltA, and OR731328 to OR731333 for groEL.

In order to investigate the distribution of A. bovis in China, we collected the collection location and host information of all A. bovis sequences in the NCBI nucleotide database. Anaplasma bovis was found to be distributed in most provinces, municipalities and autonomous regions throughout China (Figure 3(A)). Anaplasma bovis had been detected in 19 mammalian hosts and 12 tick species (Table S4), the A. bovis sequences of which were used in the downstream phylogenetic analysis. In the phylogenetic trees based on partial 16S rRNA, gltA and groEL sequences, A. bovis was always divided into two clusters (Figure 3(B–D)). Anaplasma. bovis strains involved in this study were located in different clusters, and even the 16S rRNA, gltA and groEL sequences from the same specimen were located in different clusters (Figure S3), indicating that A. bovis had high genetic diversity.

Figure 3. Records of A. bovis from NCBI nucleotide database distribution and phylogenetic analysis. (A) China map depicted A. bovis distribution situation. Colour of circle represent host of this record. Red rhombi represent records from this study. (B) Phylogenetic tree based on 16S rRNA. Colour of circle represent host of this gene. Red rhombi represent records from this study (C) Phylogenetic tree based on gltA gene. (D) Phylogenetic tree based on groEL gene.

Discussion

As A. bovis is an intracellular bacterium, it is difficult to isolate and culture. Therefore, we used the next-generation sequencing technology to directly sequence from the goat blood positive for A. bovis, and obtained the whole genome sequences of A. bovis with a completeness of 95.30%. Of all Anaplasma, A. bovis had the smallest genome. Phylogenetic analysis based on the whole genome revealed that A. bovis was closely related to A. platys and A. phagocytophilum, while it is far related to erythrocytic Anaplasma. Compared to other Anaplasma species, A. bovis lacked 51 orthogroups, but had 116 unique orthogroups. The virulence factors of A. bovis were significantly less than those of A. phagocytophilum, suggesting less pathogenicity of A. bovis. The infection rate and genetic diversity of A. bovis in Shandong goats were high, which is worthy of further investigation.

In terms of basic genome information, the genome size of A. bovis is 1.05 Mbp and the GC content is 42.94%, which is the smallest genome in Anaplasma and relatively low GC content in Anaplasma. Compared with free-living bacteria, a homogenous host environment of intracellular bacteria may create nutritional limitations that require less energy utilization. The synthesis of GTP and CTP is energy-intensive, so bacterial symbionts have genomes which were rich in A + T [41]. For intracellular bacteria, the early stage of adaptation to the intracellular lifestyle is accompanied by the loss of a large number of genes, the production of pseudogenes, and the reduction of genome size [42].

In protein annotation, proteins of A. bovis in the category cell wall/membrane/envelope biogenesis were obviously lower than other species. Proteins in cell wall/membrane/envelope biogenesis usually play a role in bacteria survival in extreme environments. For example, relA-dependent (p)ppGpp-mediated stringent response activates when bacteria face nutritional starvation. This can up-regulate genes of translation cell wall/membrane/envelope biogenesis [43]. Compared with Vibrio isolated from epipelagic, hadopelagic strains encountered high pressure, low temperature, and high inorganic salt content. It possesses more genes in cell wall/membrane/envelope biogenesis [44]. Desiccation resistance genes of Salmonella Typhimurium under desiccation stress related to cell wall/membrane/envelope biogenesis [45]. Thus, our findings suggest that A. bovis may be less able to cope with extreme environments than other Anaplasma species.

Virulence factors play an important role in bacterial establishment of infection and disease processes [46]. The virulence factors of Anaplasma can be divided into three groups: invasion, effector delivery system and immune modulation [47]. The surface exposure of VirB6-3 and VirB6-4 is associated with the growth of A. phagocytophilum in human and tick cells [48]. AnkA is an effector protein of T4SS that can bind a large number of DNA and proteins and may be involved in regulating neutrophils survival and replication. Phosphorylated AnkA interacts with SHP-1 to interfere with host cell signalling and may promote the survival of A. phagocytophilum [49–51]. OmpA facilitates pathogen invasion of mammalian and tick host cells, and anti-OmpA serum treated A. phagocytophilum reduced the ability to infect HL-60 cells [52,53]. In mammalian cells, Ats-1 can promote infection by inhibiting apoptosis and delivering autophagic cargo that can be used as nutrients by pathogens [49,54,55]. In our study, VirB6, AnkA, OmpA, ats-1 were very different from the others. This may indicate that the pathogenic effects of A. bovis may be affected in some way. Although the results of virulence factor comparisons suggest that A. bovis may be less pathogenic, its wide range of animal hosts and vectors, as well as the emergence of human infections in recent years, suggest that its public health and veterinary significance cannot be ignored.

In this study, the infection rate of A. bovis among goats in Shandong Province, China was nearly 80%, which was higher than those in Shaanxi, Zhejiang and Sichuan [7,56,57]. However, considering the small sample size of this study, the actual infection rate needs to be further investigated with large samples to clarify the significance of goats as the host of A. bovis. It is interesting that the epidemic strains of A. bovis in different regions of China are obviously divided into two clusters. In this study, the A. bovis from the same region has a higher genetic diversity. The higher prevalence of A. bovis among domestic animals may be due to its milder symptoms, which make it difficult to detect diseased animals, allowing them to survive and spread A. bovis through tick vectors. These findings suggest that field surveys and surveillance should be enhanced to prevent the transmission of the tick-borne pathogen among animals and potentially to humans [14].

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study was supported by the Natural Science Foundation of China (82103897, L.Z.; 81621005, W.-C.C.), National Key Research and Development Program of China (2021YFC2302001, X.-M.C.), the Natural Science Foundation of Shandong Province, China (ZR2020QH299, L.Z.).