Bartonella Prevalence and Genome Sequences in Rodents in Some Regions of Xinjiang, China

Rodent-borne Bartonella species have been associated with zoonotic diseases. Bartonella species such as Bartonella elizabethae, Bartonella grahamii, and Bartonella tribocorum can cause disease in humans. Humans can be infected by blood-sucking arthropods through the scratches and bites of an infected reservoir host or via contact with infectious rodents. ABSTRACT In this study, we investigated Bartonella infection and its genetic diversity in rodents in Beitun, Xinjiang Uygur Autonomous Region, China. Small mammals were captured using snap traps at four sampling sites in 2018. Spleen and liver tissues were collected and cultured to isolate Bartonella strains. Whole-genome sequencing was performed on the strains identified as Bartonella by gltA gene PCR, and the average nucleotide identity (ANI) of the genomes was calculated by using FastANI v1.33. Phylogenetic trees were constructed for the samples positive for Bartonella spp. by the gltA PCR assay based on 1,290-bp gltA genes, 2,903-bp rpoB genes, and core-genome single nucleotide polymorphisms (SNPs). Among 66 rodents, 11 were positive for Bartonella, with an infection rate of 16.67%. The rodent infection rates in different tissues (χ2 = 2.133; P = 0.242), species (χ2 = 9.631; P = 0.141), and habitats (χ2 = 4.309; P = 0.312) did not show statistical differences. Bartonella spp. isolated from the rodents were phylogenetically divided into six clades (two different Bartonella species were detected in two rodents). By comparing phylogenetic trees based on gltA genes, rpoB genes, and SNPs, we found that the topological structures of several evolutionary trees are different. However, the Bartonella strains isolated in this study were clustered into six clusters in different phylogenetic trees. Broad distributions and high genetic diversity of Bartonella strains were observed among rodents in Beitun, Xinjiang. IMPORTANCE Rodent-borne Bartonella species have been associated with zoonotic diseases. Bartonella species such as Bartonella elizabethae, Bartonella grahamii, and Bartonella tribocorum can cause disease in humans. Humans can be infected by blood-sucking arthropods through the scratches and bites of an infected reservoir host or via contact with infectious rodents. Xinjiang is one of the provinces with the most abundant species of Bartonella in China, but there are few reports about the prevalence of Bartonella in the Beitun area. This research aims to investigate the occurrence and prevalence of Bartonella infection in rodents at these sampling sites and provide a basis for the prevention and control of rodent Bartonella species in Beitun and the surrounding areas of Xinjiang.

To further verify the relationship of the Bartonella isolates from Beitun rodents with recognized species of Bartonella, we performed phylogenetic analyses based on the full-length gltA and rpoB genes. In the phylogenetic trees of the sequences (Fig. 2), 13 Bartonella sequences were distributed into six clades. The gltA and rpoB sequences obtained in this study were all clustered into the clade of Bartonella species associated with rodents. AU4XJBT, AU18XJBT, and CM31XJBT were placed into the same phylogenetic group as B. grahamii as4aup. However, in the phylogenetic tree constructed by using the 2,903-bp rpoB gene, the bootstrap value between AU4XJBT, AU18XJBT, and CM31XJBT was 65%, which could not firmly indicate that the three strains belonged to the same species (Fig.  2C). Therefore, AU4XJBT and the latter two strains were artificially classified into cluster I and cluster II, respectively. Cluster III consisted of AU5XJBT, AU15XJBT, and AU16XJBT. Cluster IV was composed of AU55XJBT, which is closest to B. elizabethae. ML70XJBT and MM73XJBT.G were classified as belonging to cluster V. The four remaining strains (ML69XJBT, ML70XJBT.G, ML71XJBT, and MM73XJBT) were divided into cluster VI. By comparing the topological structures of the three evolutionary trees, the Bartonella isolates obtained in this study were divided into six clusters, but there were differences in the details of the evolutionary trees ( Fig. 2A to C). To clarify the Bartonella species to which the 13 strains belonged, we studied the genomes of the Bartonella strains obtained in this study. Genomic features of the genus Bartonella. The genomic features of these strains, including the G1C content, genome size, and number of coding sequences (CDSs), are presented in Table 3. Briefly, the G1C content of the genomes of these strains ranged from 37.68% to 38.93%, with an average of 38.24%. The average genome size of the 13 Bartonella strains in this study was 2.50 Mb (range, 2.14 to 2.89 Mb), with numbers of CDSs ranging from 1,998 to 2,347. All of the strains obtained in this study contained one predicted tmRNA (transfer-messenger RNA), but the predicted number of tRNAs was between 39 and 44, and the predicted number of rRNAs was between 3 and 4. The number of rRNAs predicted for CM31XJBT alone was 4 (Table 3). Under a threshold of 85%, the numbers of core genes and variable genes identified by pangenomic analysis of 52 strains (13 in this study and 39 in previous publications) were 129 and 18,446, respectively. Core-genome single-nucleotide polymorphism (SNP) analysis was carried out using the snp-sites command. A total of 20,102 SNP sites were selected from the core genome, and these sites were integrated for the construction of the SNP evolutionary tree.
Genomic analysis. The average nucleotide identity (ANI) values of the 13 strains obtained in this study and 39 strains downloaded from the National Center for Biotechnology Information (NCBI) database for the genus Bartonella were used to assess the overall genome similarity ( Fig. 3; see also Table S1 in the supplemental material). The current standard for a    (Table S1). We used phylogenomic analysis to examine the evolutionary relatedness among the 13 Bartonella strains in this study and 39 other Bartonella strains (Fig. 4). The genomic assembly quality control of the strains used to construct the phylogenetic tree met the requirements (genomes with N 75 values of $10,000 bp and #500 undetermined bases per 100,000 bases) (14). Fewer SNPs were observed for core genes (129) than for total genes (18,575). The phylogenetic tree based on SNPs showed that AU4XJBT was clustered with B. grahamii. Three strains isolated from A. uralensis (AU5XJBT, AU15XJBT, and AU16XJBT) were clustered. AU55XJBT from A. uralensis had the highest similarity to B. elizabethae and B. kosoyi. The B. krasnovii clade included ML70XJBT from M. tamariscinus and MM73XJBT.G from M. meridianus. This analysis, based on SNPs in genomic sequences, showed that the ML69XJBT, ML70XJBT.G, ML71XJBT, and MM73XJBT isolates from Meriones were grouped (Fig. 4). We observed 13 Bartonella strains in this study, which were divided into six clusters using core-genome phylogenetic analysis; this was similar to the topology of the phylogenetic trees constructed with the gltA genes and rpoB genes.

DISCUSSION
In Beitun and surrounding areas, six species of Bartonella were detected in four rodent species. Among them, B. grahamii was reported to be able to cause neuroretinitis in humans  Table S1 in the supplemental material.

Prevalence and Genome Sequences of Bartonella
Applied and Environmental Microbiology (15) and bilateral retinal artery branch occlusions (16). B. tribocorum has been reported in fleas of domestic dogs (17) and was implicated as the causative agent of human infections (18). In this study, the prevalence of Bartonella species in infected rodents varied from 0 to 50.0%, and the highest rate was detected in M. meridianus. Two or more Bartonella species were detected in A. uralensis. At least two species of Bartonella were isolated from both M. tamariscinus and M. meridianus. A B. grahamii-like isolate was detected in both the liver and spleen of a gray hamster. There was no statistically significant differences in the prevalences of Bartonella infection among the various rodents. According to the sequencing of the 320bp gltA gene fragment, two strains from A. uralensis were more than 97% identical to B. grahamii, and one strain from C. migratorius was 99.69% similar to B. grahamii (Table 2). These results suggested that there might be intraspecies variations in B. grahamii from different host species; a similar phenomenon was previously reported on Heixiazi Island in northeastern China (19). In this study, we found that rodents 70 and 73 were coinfected by two Bartonella species using culture experiments. This suggests that we need to pay attention to this phenomenon in follow-up pathogen surveillance to prevent and control Bartonella disease better.
In this study, the species of Bartonella isolated were confirmed at both the gene and genome levels. Both gene-based phylogenetic analysis and SNP-based phylogenetic trees showed that the 13 Bartonella strains in this study were divided into 6 species ( Fig.  2 and 4). However, in different phylogenetic trees, the closest known Bartonella species (in the NCBI database) of the 13 strains isolated in this study were different, especially those with ANI values of ,95%. In general, an ANI value of ,95% is used as a cutoff point for determining new species (11,12,20). In the ANI analysis of the genomic sequences of the Bartonella strains, all strains isolated from Xinjiang had ANI values of ,95% compared to the Bartonella sequences in GenBank, except for the AU4XJBT strain. A previous study from Brazil reported that most Bartonella species showed ANI values of between 75% and 90%, and a cutoff point of 95% , ANI , 97% for the classification of new Bartonella subspecies was proposed (21). According to the existing ANI threshold, AU18XJBT and CM31XJBT do not belong to the same species. Similarly, ML70XJBT.G and the other three strains in cluster VI are also not the same species. Some inconsistencies were observed between these results and the phylogenetic trees based on the gltA gene, the rpoB gene, and the core-genome SNPs. On the one hand, the reasons might be the diversity and complexity of the Bartonella genome (18,446 variable genes). On the other hand, it may be due to gaps in the assembled genome. More Bartonella sequences are needed to verify whether there are different results of the phylogenetic analyses of the draft genome and the complete genome.
In this study, we could not accurately identify four strains of Bartonella (ML69XJBT, ML70XJBT.G, ML71XJBT, and MM73XJBT). They were highly similar to "Candidatus Bartonella gerbillinarum" (99.69% to 100%) by 320-bp gltA gene alignment. Comparison of the fulllength gltA genes (1,296 bp) showed that B. alsatica had the highest similarity (90.28% to 90.51%) with these strains. B. taylorii had the highest similarity (91.04% to 91.26%) with these strains by comparison of the full-length rpoB genes (4,152 bp). All values were below the threshold of 97%. Phylogenetic trees based on SNPs in the core genome showed that the four strains mentioned above were closest to B. florencae R4. At the whole-genome level, the four strains had the highest ANIs (85.95% to 86.12%) with B. florencae R4, which did not reach the cutoff value that could accurately identify the species. These results indicate that ML69XJBT, ML70XJBT.G, ML71XJBT, and MM73XJBT may potentially be new species of Bartonella. The ANI values of ML70XBT.G and the other three strains were between 93.65% and 93.76%, suggesting that ML70XBT.G might be a subspecies.
In summary, based on the isolation results, the livers and spleens of infected animals had similar Bartonella burdens. These and previous observations suggest that the long-term persistence of Bartonella in these tissues may damage these organs (22). Comparative core-genome SNP, ANI, and gene analyses confirm the taxonomic placement of the Bartonella strains isolated in this study with other Bartonella spp. This study demonstrated the existence of multiple Bartonella infections in rodents in Beitun, Xinjiang. At the same time, Bartonella coinfection was found in rodents in this region. We suggest that more attention should be paid to coinfection in rodent-associated Bartonella surveillance in order to better understand the prevalence and genetic diversity of Bartonella species in rodents.

MATERIALS AND METHODS
Specimen collection and identification of rodents. Rodents were captured using the trap-at-night method in June 2018 in Beitun and surrounding areas of the Xinjiang Uygur Autonomous Region (longitude, 87°539 to 87°979; latitude, 47°279 to 47°579). Traps baited with peanuts were set up in three regions consisting of five habitats, an Elaeagnus angustifolia plantation, a shelterbelt forest, tamarisk trees, grassland, and a semidesert. Sampling regions included Altay City, Beitun City, and Fuhai County, including four sampling sites (Fig. 1). Morphological and DNA barcoding techniques were used to identify the species of the trapped animals. Animals were accurately identified by the amplification of the mitochondrial cytochrome c oxidase subunit I (COI) gene with the universal primers shown in Table 4 (23). The collection time, site, habitat, species, gender, weight, head-body length, and tail length were recorded. The liver and spleen were collected into 1.5-mL Nunc CryoTubes and stored in liquid nitrogen prior to processing and laboratory testing.
Isolation and purification of Bartonella organisms. Approximately 10 mg of organ samples (liver or spleen) was added to 2-mL centrifuge tubes with 100 mL of tryptic soy broth (TSB) (BD) and some 1.2mm grinding beads and disrupted using a Dounce-Potter homogenizer under sterile conditions, and the suspensions were inoculated onto tryptic soy agar (TSA) (BD) plates containing 5% defibrinated sheep blood. Plates were kept in a humid environment with 5% CO 2 at 36.5°C for 20 days. Next, we selected suspected Bartonella organisms for purification and passage 2 to 4 times.
Material from a single colony was suspended in 100 mL of deionized water, heated for 30 min at 100°C, and cleared by centrifugation for 5 min at 6,000 Â g at 4°C. The supernatant was used as a DNA template for PCR. The initial screening of Bartonella DNA was performed using PCR to detect a portion of the citrate synthase gene (gltA) with the primers BhCS.781p and BhCS.1137n (Table 4), as described previously (24). The pure bacterial cultures were collected in brain heart infusion (BHI) broth (BD) supplemented with 30% glycerol at 270°C (19,25).
Genome sequencing, assembly, and annotation. To further characterize the positive samples in this study, genomic DNA from positive samples was extracted using the Wizard genomic DNA purification kit (Promega, Madison, WI, USA) according to the manufacturer's instructions. Genomic DNA was extracted by using SDS (26), followed by agarose gel electrophoresis to determine the purity and integrity of the DNA, and quantification was performed by using a Qubit 2.0 fluorometer. DNA samples quantified by electrophoresis were randomly interrupted by using a Covaris ultrasonic breaker at a length of about 350 bp. After the processing of complete DNA fragments, the entire library preparation was completed using the NEBNext Ultra DNA library prep kit for Illumina (New England BioLabs, USA). The Qubit 2.0 instrument was used for the initial quantification of the library, the library was diluted to 2 ng/mL, and the inserts in the library were subsequently detected using an Agilent 2100 instrument. Finally, the effective concentration of the library was accurately quantified by using a quantitative PCR (qPCR) method to ensure the quality of the library. After library inspection, different libraries were sequenced on the Illumina NovaSeq PE150 (150-bp paired-end) platform according to the effective concentration and the target amount of data. All sequencing depths were about 100-fold. Readfq (version 10) was used to filter the raw data and obtain effective data (clean data) for downstream analysis. The specific processing steps were as follows: (i) reads with .40% low-quality bases (mass value of #20) were filtered out, (ii) reads containing .10% N bases were filtered out, (iii) reads with overlaps and adapter sequences exceeding a certain threshold (the default was 15 bp) were removed, and (iv) host-derived reads and duplications were filtered out. All goodquality paired-end reads were then de novo assembled using SPAdes v3.8.0 (27). The sequenced genomes were annotated using Prokka v1.14.5 (28).
Phylogenetic analysis of the gltA and rpoB genes. The DNAStar Lasergene SeqMan (plug-in in the DNAStar installation package) procedure was applied for sequence splicing, and a BLAST search was performed (http://blast.ncbi.nlm.nih.gov/Blast.cgi) to determine splicing sequence similarity values. Bartonella isolates were considered similar to validated Bartonella species if the gltA fragments had $97% sequence similarity (29). The full-length gltA genes and rpoB genes of the Bartonella strains cultured in this study were screened from whole-genome sequences using localized BLAST analysis (ncbiblast-2.13.01). Gene sequences were aligned using ClustalW. A phylogenetic tree was created by using the maximum likelihood (ML) method in Mega_X_10.1.8 with 1,000 bootstrap replicates. The best evolutionary model was selected by the Models program in Mega_X_10.1.8. Reverse ACTTCTGGGTGTCCAAAGAATCA Phylogenomic analyses. Pairwise average nucleotide identity (ANI) values were estimated using FastANI v1.33 (12). The results were visualized and plotted with the R package pheatmap. The assembled Bartonella isolate genomes (this study) were compared with 39 published Bartonella genomes with valid names retrieved from the National Center for Biotechnology Information (NCBI) RefSeq database (30), but some Bartonella genomes were not included in the RefSeq database (Tables 2 and 5). At the same time, the core genomes and pangenomes of the Bartonella genomes mentioned above were compared. Single nucleotide polymorphisms (SNPs) in the core genomes were found using Roary v1.13.0 (31) analysis of these genomes. The minimum percent identity by BLASTp for Roary software was set to 85%. These SNPs were aligned with MAFFT (32) and were later used for phylogenetic analyses. The phylogenetic tree was inferred with IQ-TREE v2.1.2 by using the ML method. The best model was selected by the ModelFinder program. The number of bootstrap replicates was 1,000.
Coinfection certification experiment. When the gltA gene was used for the preliminary identification of Bartonella species obtained in this study, an interesting phenomenon was found: different species of Bartonella were cultured from the spleen and liver tissues of animals 70 and 73. To verify the presence of two Bartonella coinfections in rodents 70 and 73, tissue samples from these two animals were recultured. In this Bartonella isolation experiment, we selected as many monoclonal bacteria as possible from the primary culture medium for passage culture. Next, the appropriate numbers of colonies were selected from the purified culture medium to prepare a DNA template for PCR, and the gltA gene was then amplified for identification.
Statistical analysis. R v4.1.2 was used for statistical analysis, with the significance level set at a P value of 0.05.
Downloading of publicly available assemblies. The assembly genome sequences of validly named Bartonella species were downloaded from the NCBI public database on 1 July 2022. All publicly available assemblies were subjected to quality control by using Quast v5.0.2 (33). Genomes with N 75 values of ,10,000 bp and .500 undetermined bases per 100,000 bases were discarded (14). Table 5 shows the Bartonella species downloaded from the NCBI that were included in this genomic analysis. Data availability. The novel genome sequences of the Bartonella strains have been submitted to the NCBI GenBank database under accession numbers JAQNDQ000000000, JAQITG000000000, JAQITF000000000, JAQNDP000000000, JAQITE000000000, JAQITC000000000, JAQITD000000000, JAQITB000000000, JAQISX000000000, JAQITA000000000, JAQISZ000000000, JAQISY000000000, and JAQISW000000000.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, XLSX file, 0.02 MB.