The spleen bacteriome of wild rodents and shrews from Marigat, Baringo County, Kenya

Sample acquisition

The animal use protocol for this study was reviewed and approved by Walter Reed Army Institute of Research under protocol number AP-12-001, KEMRI IACUC #2208, and National Museums of Kenya NMK/SCom2013/08. Archived spleen tissues from a previous study were used (Masakhwe et al., 2018). Briefly, the rodent survey was conducted between 8th and 14th of June 2017. Each night, Sherman collapsible rodent traps were set up in a variety of habitats including cropped fields, grass fields, gardens, orchards, and around buildings. Trapped animals were removed and after euthanasia, the surface of the rodent abdomen was cleaned with 70% ethanol before necropsy. About 100 mg of the spleen samples were collected, stored in screw cap plastic vials before placing them in shippers containing liquid nitrogen. The tissue biopsies were later transferred to the laboratory where they were stored at −80 °C.

Genomic DNA isolation from the spleen tissues

Total genomic DNA was isolated from ∼10 grams of spleen tissues using QIAamp DNA Mini Kit (Qiagen, Valencia, CA, USA) as recommended by the manufacturer. DNA was eluted in 100 µL and stored in a −80 °C freezer until required.

Molecular taxonomy to aid identification of individual animal species

The primers used to amplify the cytb gene are listed in Table S1. Where amplification could not be achieved with the primary PCR primers, a secondary nested PCR was performed using internal primers L14749 and H14896 as described by (Nicolas et al., 2012). PCR was performed in a 25 µL reaction volume containing 2 µL of DNA template, 0.2 µM of each primer, 0.05 U of MyTaq polymerase and 5 µL of 5X MyTaq buffer (Bioline, UK). Amplification was performed in an Eppendorf Master cycler pro 384 (Eppendorf, USA) with an initial denaturation step of 94 °C for 2 min, 30 cycles of 94 °C for 30 s, 52 °C for 30 s, and 72 °C for 1 min followed by a final extension step at 72 °C for 10 min. Amplicons were visualized on a 2% agarose gel (Thermo Fisher Scientific, CA, USA) stained with GelRed (Biotium, Australia). PCR products from positive samples were purified using AMPure XP beads (Beckman Coulter, CA, USA). The forward and reverse primers were used for cycle sequencing by the Big Dye Terminator Cycle Sequencing Kit v 3.1 (Applied Biosystems, CA, USA). Reaction products were then purified using CleanSEQ beads (Beckman Coulter, CA, USA) and sequenced on a 3130 Genetic Analyzer (Applied Biosystems CA, USA).

Spleen bacterial diversity by 16S rRNA deep sequencing

The bacteria DNA in the spleen biopsy was amplified with primers that target the 16S rRNA V3–V4 hyper variable region (Table S1). The primers were tagged with Illumina sequencing adapters as previously described (Klindworth et al., 2013). A non-target control sample (PCR water) was used to track laboratory contaminants. Amplification was carried out in a 25 µL reaction consisting of 12.5 µL of 2X NEB Next PCR master mix (New England Biolabs, MA, USA), 0.2 µM of each primer and 2.5 µL of DNA template. Amplification was performed on the Eppendorf Master cycler pro 384 (Eppendorf, Humberg, Germany) with an initial denaturation at 95 °C for 3 min followed by 25 cycles of 95 °C for 30 s, 55 °C for 30 s, 72 °C for 30 s and a final extension of 72 °C for 5 min. Amplicons were thereafter purified using AMPure XP beads according to the manufacturer’s instructions (Beckman Coulter Genomics, Brea, CA, USA).

A dual indexing PCR to allow multiplexing of samples was done in a 50 µL reaction consisting of 5  µL of purified amplicons, 5  µL of each Nextera XT i7 and i5 Index Primer (Illumina, USA), 25 µL of NEBNext High-Fidelity 2X PCR Master Mix (New Englands BioLabs, MA, US) and 10 µL of PCR grade water (Thermo Fisher Scientific, Waltham, MA, USA), with thermocycling at 95 °C for 3 min, followed by 12 cycles of 95 °C for 30 sec, 55 °C for 30 sec, and 72 °C for 30 sec, and a final extension at 72 °C for 5 min. The indexed amplicon libraries were purified with AMPure XP beads according to the manufacturer’s instructions (Beckman Coulter Genomics, Brea, CA, USA), and then quantified on Qubit Fluorometer 2.0 using Qubit dsDNA HS assay kit according to the manufacturer’s protocol (Thermo Fisher Scientific, Waltham, MA, USA). Libraries were normalized to a concentration of 4 nM and then pooled. The pooled samples were denatured and diluted to a final concentration of 12 picomoles, then spiked with 5 % PhiX (Illumina, USA) as a sequencing control. Samples were sequenced on MiSeq platform (Illumina, USA) using MiSeq 600 cycle reagent kit V3 (Illumina, USA).

Data analysis

Forward and reverse nucleotide sequences from the cytb gene were quality checked and assembled into contigs using CLC main work bench version 7.5 (CLC Inc, Aarhus, Denmark). Identification of species within the orders Rodentia and Eulipotyphyla was done by querying the assembled contigs against the nucleotide database (GenBank) using the Nucleotide Basic Local Alignment Search Tool (BLASTn) (Altschul et al., 1990).

To determine the phylogenetic placements of the small mammals, reference sequences corresponding to the four small mammal species i.e., Acomys wilsoni, Mastomys spp, Crocidura spp and Rattus rattus were downloaded from GenBank (http://www.ncbi.nlm.nih.gov/genbank/). Multiple sequence alignments were performed in MEGA version 7 (Kumar, Stecher & Tamura, 2016) and the General Time Reversal model with gamma distribution and invariable sites (GTR+ Γ+I) was determined as best model for the analysis. The alignments were then used to infer phylogenetic relationship using the maximum likelihood treeing method. Branch support was calculated as 1,000 bootstrap replicates.

For the 16S rRNA deep sequencing, the Miseq sequences output were de-multiplexed and adapters trimmed using the Miseq reporter software version 2.6.3 (Illumina, USA). Mothur version 1.35 pipeline was used for paired end reads contig assembly, sequence quality filtering, chimera removal and taxonomic identification (Schloss et al., 2009). In brief, contigs containing ambiguous bases, and those with lengths ≤389 bp and ≥500 bp were discarded. The SILVA database was customized for the V3-V4 region and used for sequence alignment, followed by merging sequences that were not more than 2 bp different from each other using the pre-cluster command in Mothur (Quast et al., 2013). Unassigned OTUs and those assigned to Chloroplast, Mitochondria, Archaea, and Eukaryota were discarded. Samples with less than 1,000 (n = 10) sequences were discarded as small library sizes are often a confounding factor that obscures biologically meaningful results (Weiss et al., 2017). Taxa contained in the non-template control were censored from sample dataset as proposed by Davis et al. (2018).

Statistical analysis and visualization were performed using R software environment (version 4.1, R Core Team, 2016), with the following add-on statistical packages: ‘phyloseq’ (McMurdie & Holmes, 2013), ‘vegan’ and ‘ggplot2’ (Oksanen et al., 2014; Wickham, 2009). Rarefaction was performed to down sample the data for alpha diversity analysis using the rarefy_even_depth command in phyloseq with replacement. This was done to account for unequal sequencing between samples (Weiss et al., 2015). The rarefied data was used to compute the Shannon diversity index (Shannon, 1948), by first determining the bacterial diversity for each animal, and then aggregated to genus level before computing the mean.

Small mammal speciation by cytb gene

A total of 54 spleen tissues from wild caught small mammals were used for molecular taxonomy and microbiome profiling. The small mammals had previously been classified based on morphological characteristics into two major orders, Rodentia and Eulipotyphyla (Masakhwe et al., 2018). The order Rodentia (38/54 =70.4%) had the following species: Acomys wilsoni (n = 21), Lophuromys sikapusi (n = 1), Mastomys natalensis (n = 8), Rattus rattus (n = 2), Arvicanthis niloticus (n = 6). The order Eulipotyphyla (16/54 =29.6%) had Crocidura olivieri (n = 2), and other un-speciatable Crocidura spp (n = 14).

Using cytb gene, 49/54 samples generated usable sequences. Homology searches with BLASTn against the GenBank database classified the small mammals into two main orders: Rodentia (77.6%) comprising Acomys (21/49), Mastomys (8/49), Arvicanthis (6/49) and Rattus (3/49), and Eulipotyphyla (22.4%), comprising only Crocidura (11/49).

Of the 49 samples sequenced, only 46 could be used for phylogenetic analysis. These were compared against 20 representative species in a maximum likelihood tree. As shown in Fig. 1, the shrews (11/46) clustered with Crocidura somalica in a well-supported clade (99%). The rodent species were more diverse. 20/46 samples branched with Acomys wilsoni (100%), 3/46 samples branched with Rattus rattus, 6/46 branched with Arvicanthis niloticus and 6/46 branched with Mastomys natalensis.

Diversity of bacteria in the small mammals

Sequences from the 54 spleen samples yielded 8,621,961 raw sequence contigs. After quality filtering, collapsing duplicate sequences, removing chimeras and non-bacterial sequences, 2,778,822 sequences were considered suitable for further analysis. On querying the SILVA rRNA database, the sequences were grouped into 182 taxa at 97% sequence similarity. Of the 54 spleen samples examined, 10 that had small library sizes (<1000 sequences) were dropped from downstream analysis. The remaining 44 samples had a total of 2,778,038 sequences. After contamination screening a total of 7 OTUs were identified as contaminants by the prevalence method in the “decontam” R package (see Table S2 for list of contaminants) and were removed from the spleen sample dataset.

Figure 2 shows the 17 bacteria phyla that were detected. Proteobacteria was the most abundant and contributed 35,122 of 51,731 total contigs (67.9%). Other phyla included Actinobacteria at 16.5%%, Firmicutes (5.5%), Chlamydiae (3.8%), Chloroflexi (2.6%) and Bacteroidetes (1.3%). Less abundant phyla (<1%) included Acidobacteria, Verrucomicroba, Planctomycetes, Fusobacteria, Deinococcus-Thermus, Armatimonadetes, Gemmatimonadetes, TM7, OD1, Spirochaetes and SR1.

Taxonomic assignment at phylum and genus level for the 182 OTUs is shown in Fig. 3. Of these genera, Bartonella was the most abundant with sequence yield of 45.6% (19,872 out of 43,622 total contigs). Because some of the bacteria could not be classified beyond family level, the total contigs at genera level were fewer (43,622) than those at phylum level (51,731 contigs). Compared to Bartonella, the relative abundance of other bacteria in the Proteobacteria was low: Anaplasma (8.0%), Delftia (3.8%), Methylobacterium (3.5%), Coxiella 2.6%, Bradyrhizobium (1.6%) and Achromobacter (1.2%), Acinetobacter (1.1%). Potentially pathogenic bacteria from other phyla occurred at <1% and included Ehrlichia, Rickettsia and Brucella.

Figure 4 shows a heat map of 44 spleen samples from the different small mammal species clustered against the bacteria genera with abundance >5%. Of the 44, 34 had Bartonella sequences. Of these, 20 (7 from Mastomys, 5 Crucidura, 4 Arvicanthus, 3 Acomys and 1 Ratus) contributed to the over-representation of Bartonella sequences in the heat map. Interestingly, the spleen samples with over-representation of Bartonella tended to have fewer or no other bacteria genera.

The Shannon diversity index ranged from 0.3204 to 3.5135, indicating a broad variation in the bacterial diversity between samples. The mean bacteria diversity in Acomys was 2.86 compared to 1.77 for Crocidura, 1.44 for Rattus, 1.40 for Arvicathis and 0.60 for Mastomys (Fig. 5). Using one-way ANOVA followed by the Tukey HSD (Honest Significant Difference) posttest, the bacteria diversity was significantly higher in Acomys compared to Mastomys (p = 0.016).

Disclaimer

Material has been reviewed by the Walter Reed Army Institute of Research. There is no objection to its publication. The opinions or assertions contained herein are the private views of the authors, and they are not to be construed as official, or as reflecting true views of the Department of the Army or the Department of Defense.