-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathThai_bats_supplement.Rmd
208 lines (173 loc) · 20.1 KB
/
Thai_bats_supplement.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
---
title: "Supplementary Material"
subtitle: Diversity and phylogenetic relationships among Bartonella strains from Thai bats
author:
- affiliation: Division of Vector-Borne Diseases, Centers for Disease Control and
Prevention, Fort Collins, CO 80521, USA; Department of Biology, Colorado State
University, Fort Collins, CO 80523, USA
name: Clifton D. McKee
- affiliation: Division of Vector-Borne Diseases, Centers for Disease Control and
Prevention, Fort Collins, CO 80521, USA
name: Michael Y. Kosoy, Ying Bai, Lynn M. Osikowicz
- affiliation: Division of High-Consequence Pathogens and Pathology, Centers for Disease
Control and Prevention, Atlanta, GA 30333, USA; National Wildlife Research Center,
USDA/APHIS/Wildlife Services, Fort Collins, CO 80521, USA
name: Amy T. Gilbert
- affiliation: Division of High-Consequence Pathogens and Pathology, Centers for Disease
Control and Prevention, Atlanta, GA 30333, USA
name: Richard Franka
- affiliation: Faculty Sciences and Public Health, Rajapruk University, Nonthaburi
11130, Thailand
name: Sumalee Boonmar
- affiliation: LYSSA LLC, Lawrenceville, GA 30044, USA
name: Charles E. Rupprecht
- affiliation: Center for Global Health, Centers for Disease Control and Prevention,
Atlanta, GA 30333, USA
name: Leonard F. Peruski
output:
html_document:
fig_height: 8.5
fig_width: 11
pdf_document: default
word_document: default
bibliography: "Thai_bats_supplement.bib"
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
library(reshape2)
library(mapdata)
library(maps)
library(maptools)
library(rgeos)
library(sp)
library(binom)
library(knitr)
```
#### Supplementary Materials and Methods
##### Sampling locations and bat capture
Bats were captured from locations in four provinces (Fig A): Wat Tham Phra Cave in Mueang Chiang Rai District, Chiang Rai Province; Sia Cave in Khlong Lan National Park, Khlong Lan District, Kamphaeng Phet Province; a cave in Phu Pha Man District, Khon Kaen Province; and Khao Chakan Cave, Mueang Sa Kaeo District, Sa Kaeo Province.
Bat capture followed the approach used previously [@Kuzmin2008; @Kosoy2010b; @Bai2012a]. Briefly, bat sampling sites in each province were chosen based on available information about bat roosts and observations of bat activity in the area. Bats were collected using hand nets or manually in caves while mist nets were used to capture foraging bats. Captured bats were anesthetized by intramuscular injection of ketamine hydrochloride (0.05-0.1 mg/g body weight) and euthanized under sedation in accordance with the field protocol approved by the CDC Institutional Animal Care and Use Committee; the CDC IACUC also specifically approved this study. Bats were weighed, sexed, and identified phenotypically to species using available field keys; some individuals could only be identified to the genus level. Bats were exsanguinated by cardiac puncture and blood was stored in sterile plastic tubes. Blood samples were transported from the field on dry ice.
```{r Thailand_map, echo=F, include=T, message=F, warning=F}
Thailand1=readRDS("THA_adm1.rds")
mycolors <- rep("lightgrey", dim(Thailand1@data)[1])
mycolors[c(11, 15, 17, 55)] <- "black"
plot(Thailand1, col = mycolors, border = 'darkgrey')
text(x=Thailand1@polygons[[11]]@labpt[1]+1, Thailand1@polygons[[11]]@labpt[2]+1,
labels="Chiang Rai", cex=1)
lines(x=c(Thailand1@polygons[[11]]@labpt[1], Thailand1@polygons[[11]]@labpt[1]+1),
y=c(Thailand1@polygons[[11]]@labpt[2], Thailand1@polygons[[11]]@labpt[2]+0.7))
text(x=Thailand1@polygons[[15]]@labpt[1]-3, Thailand1@polygons[[15]]@labpt[2],
labels="Kamphaeng Phet", cex=1)
lines(x=c(Thailand1@polygons[[15]]@labpt[1], Thailand1@polygons[[15]]@labpt[1]-1.2),
y=c(Thailand1@polygons[[15]]@labpt[2], Thailand1@polygons[[15]]@labpt[2]))
text(x=Thailand1@polygons[[17]]@labpt[1]+3, Thailand1@polygons[[17]]@labpt[2]+2,
labels="Khon Kaen", cex=1)
lines(x=c(Thailand1@polygons[[17]]@labpt[1], Thailand1@polygons[[17]]@labpt[1]+3),
y=c(Thailand1@polygons[[17]]@labpt[2], Thailand1@polygons[[17]]@labpt[2]+1.7))
text(x=Thailand1@polygons[[55]]@labpt[1]+1.5, Thailand1@polygons[[55]]@labpt[2]-0.5,
labels="Sa Kaeo", cex=1)
lines(x=c(Thailand1@polygons[[55]]@labpt[1], Thailand1@polygons[[55]]@labpt[1]+1.5),
y=c(Thailand1@polygons[[55]]@labpt[2], Thailand1@polygons[[55]]@labpt[2]-0.2))
```
**Fig A. Map of bat sampling locations in Thailand.** Bats were sampled in four provinces: Chiang Rai, Kamphaeng Phet, Khon Kaen, and Sa Kaeo. Counts of bat species captured in each location are shown in Figure S2.
##### *Bartonella* spp. culturing
Blood samples from bats were cultured following previously published protocols for *Bartonella* bacteria [@Kosoy1997]. Briefly, bat blood was diluted 1:4 in brain heart infusion broth with 5% fungizone (amphotericin B) and 100 µl of the sample was placed on agar plates supplemented with 10% rabbit blood. Plates were then incubated at 35 ºC with 5% $\text{CO}_{2}$ for up to five weeks, checking periodically for growth. Bacterial colonies that resembling those typical for *Bartonella* were passaged onto new plates to obtain pure cultures. Morphologically unique colonies obtained from the same sample were subpassaged in an attempt to find possible *Bartonella* coinfections. All isolates were collected in 10% glycerol solution. Cultures were then sent to the CDC Division of Vector-Borne Diseases in Fort Collins, Colorado for genetic characterization.
##### Priors used in Bayesian phylogenetic analyses
For all phylogenetic analyses using BEAST v1.8.4 on the five genetic loci analyzed (*ftsZ*, *gltA*, *nuoG*, *rpoB*, and ITS), the following default, diffuse priors were used for the GTR+Γ+I and birth-death models:
* A-C substitutions ~ gamma(0.05, 10), initial value = 1 </br>
* A-G substitutions ~ gamma(0.05, 20), initial value = 1 </br>
* A-T substitutions ~ gamma(0.05, 10), initial value = 1 </br>
* C-G substitutions ~ gamma(0.05, 10), initial value = 1 </br>
* G-T substitutions ~ gamma(0.05, 10), initial value = 1 </br>
* Base frequences ~ uniform(0, 1), initial value = 0.25 </br>
* Gamma shape parameter ~ exponential(0.5), initial value = 0.5 </br>
* Proportion of invariant sites ~ uniform(0, 1), initial value = 0.5 </br>
* Relative rates among partitions ~ uniform(0, 1E100), initial value = 1 </br>
* Birth-death birth rate ~ uniform(0, 1E5), initial value = 0.01 </br>
* Birth-death relative death rate ~ uniform(0, 1), initial value = 0.5 </br>
* Proportion of taxa sampled from birth-death tree ~ beta(1, 1), initial = 0.01
Previous sensitivity analyses have determined that these priors are sufficiently diffuse to have little influence on the posterior distributions of the parameters. The strict clock rate was fixed at one so that the branch lengths of the phylogenetic trees are scaled as substitutions per site.
#### Supplementary Results
##### Bat species distributions
The distribution of bat species among the sampled locations was highly variable (Fig B). Specifically, there was limited overlap in species among locations and most locations were dominated by one species. Stoliczka's trident bats (*Aselliscus stoliczkanus*, Hipposideridae), a fulvus roundleaf bat (*Hipposideros fulvus*, Hipposideridae), horshoe bats (*Rhinolophus* sp., Rhinolophidae), and a black-bearded tomb bat (*Taphozous melanopogon*, Emballonuridae) were captured in Wat Tham Phra Cave in Chiang Rai. Sia Cave in Kamphaeng Phet was dominated by intermediate roundleaf bats (*H. larvatus*), with a few great roundleaf bats (*H. armiger*) and a single horshoe bat (*Rhinolophus* sp.). The bat cave in Khon Kaen produced mostly *H. armiger*, but also some wrinkle-lipped free-tailed bats (*Chaerephon plicatus*, Molossidae) and a single roundleaf bat (*Hipposideros* sp.). Khao Chakan Cave in Sa Kaeo produced only *C. plicatus*. Not all samples produced adequate blood samples for culturing, so the total counts per species were slightly different between Fig B and Table A.
```{r species_sampling, echo=F, include=T, message=F, warning=F}
counts = read.csv("Thai_bats_counts.csv")
m.counts = melt(counts)
ggplot(data=m.counts, aes(x=Location, y=value, group=Species, fill=Species)) +
geom_bar(stat="identity") +
theme_bw(base_size=20) +
ylim(0, 80) +
ylab("Count") +
scale_fill_brewer(palette="Dark2")
```
```{r prop_test, echo=F, include=F, message=F, warning=F}
# Import data for prevalence
prevalence = read.csv("Thai_bats_prevalence.csv")
# Binomial confidence intervals
prev.conf = binom.confint(x=prevalence$Positive, n=prevalence$Count, methods="exact")
prev.conf[, 4:6] = round(prev.conf[, 4:6]*100, 1)
prev.conf$Province = prevalence$Province
prev.conf$Species = prevalence$Species
# Binomial test of proportions for locations
loc.test = prop.test(x=c(sum(prevalence[1:3, 4]), sum(prevalence[4:7, 4]), prevalence[8, 4],
sum(prevalence[9:11, 4])),
n=c(sum(prevalence[1:3, 3]), sum(prevalence[4:7, 3]), prevalence[8, 3],
sum(prevalence[9:11, 3])))
loc.conf = binom.confint(x=c(sum(prevalence[1:3, 4]), sum(prevalence[4:7, 4]), prevalence[8, 4],
sum(prevalence[9:11, 4])),
n=c(sum(prevalence[1:3, 3]), sum(prevalence[4:7, 3]), prevalence[8, 3],
sum(prevalence[9:11, 3])), methods="exact")
loc.conf[, 4:6] = round(loc.conf[, 4:6]*100, 1)
loc.conf$Province = c("Chiang Rai", "Kamphaeng Phet", "Sa Kaeo", "Khon Kaen")
# Total binomial confidence interval
conf = binom.confint(x=34, n=93, methods="exact")
conf[, 4:6] = round(conf[, 4:6]*100, 1)
# Binomial test of proportions for species (aggregated across locations)
s.prevalence = prevalence[order(prevalence$Species),]
sp.test = prop.test(x=c(s.prevalence[1, 4], sum(s.prevalence[2:3, 4]), sum(s.prevalence[4:5, 4]),
s.prevalence[6, 4], s.prevalence[7, 4], s.prevalence[8, 4],
sum(s.prevalence[9:10, 4]), s.prevalence[11, 4]),
n=c(s.prevalence[1, 3], sum(s.prevalence[2:3, 3]), sum(s.prevalence[4:5, 3]),
s.prevalence[6, 3], s.prevalence[7, 3], s.prevalence[8, 3],
sum(s.prevalence[9:10, 3]), s.prevalence[11, 3]))
# Binomial test of proportions for species (not aggregated across locations)
prop.test(x=prevalence[, 4], n=prevalence[, 3])
```
**Fig B. Counts of bat species captured in four provinces of Thailand.** Species were identified as close to the species level as possible. Three individuals were only identified to the genus level (*Hipposideros* and *Rhinolophus* sp.).
##### *Bartonella* prevalence and distribution of genogroups
There was some variation in *Bartonella* spp. prevalence among bat species and capture locations (Table A). In Chiang Rai, no bartonellae were cultured from *Rhinolophus* sp. or *A. stoliczkanus*, however the one *H. fulvus* individual captured was positive for *Bartonella* genogroup H3 according to its *gltA* sequence (Fig 1). In Kamphaeng Phet, *H. armiger* and *Rhinolophus* sp. showed no infection while `r prev.conf[5, 4]`% (12/23) of *H. larvatus* individuals and the single *T. melanopogon* individual were harboring bartonellae. *H. larvatus* individuals were found to be carrying *Bartonella* genogroups H1, H2, and H3 and the *T. melanopogon* was carrying *Bartonella* genogroup Tm. In Sa Kaeo, `r prev.conf[8, 4]`% (17/41) of captured *C. plicatus* individuals were carrying bartonellae; *Bartonella* genogroups Cp1, Cp2, and Cp3 were found in these bats. In Khon Kaen, *C. plicatus* and *Hipposideros* sp. individuals were negative for bartonellae, while `r prev.conf[9, 4]`% (3/14) of captured *H. armiger* individuals were infected with *Bartonella* genogroups H1 and H3. Despite the observed variation in *Bartonella* spp. prevalence, the differences among locations ($\chi^{2}$ = `r round(loc.test$statistic, 2)`, p-value = `r round(loc.test$p.value, 2)`) and species ($\chi^{2}$ = `r round(sp.test$statistic, 2)`, p-value = `r round(sp.test$p.value, 2)`) were not statistically significant, most likely due to the small sample size of bat species from each location. Larger sample sizes will be required to find significant differences among bat species and locations.
The counts of genogroups are larger than to total number of infected individuals due to the presence of multiple cultures originating from the same individual bat (Table A). The *T. melanopogon* individual (KP283) from Kamphaeng Phet produced two cultures (KP283a and KP283b), both of which were *Bartonella* genogroup Tm according to *gltA* (Fig 1); due to this similarity, only culture KP283b was chosen for multi-locus characterization. One *C. plicatus* individual (SK198) from Sa Kaeo produced two cultures (SK198a and SK198b) which were both *Bartonella* genogroup Cp3, so only culture SK198a was chosen for multi-locus characterization. One *H. armiger* individual (KK200) produced two cultures (KK200a and KK200b) which were both *Bartonella* genogroup H1, so only culture KK200a was characterized with other loci. Two *H. larvatus* individuals from Kamphaeng Phet (KP216 and KP268) produced two cultures each of genogroup H1 (KP216a, KP216b, KP268a, and KP268b), so only KP216a and KP268a was further characterized. Another *H. larvatus* bat from Kamphaeng Phet (KP287) produced multiple cultures, but only one (KP287a) was confirmed as a *Bartonella* species. Finally, *H. larvatus* KP293 from Kamphaeng Phet produced one culture (KP293a) which was genogroup H1 and a second culture (KP293b) which was genogroup H2. Only culture (KP293a) was chosen for further characterization since KP293b was so similar to culture KP277 (genogroup H2). Additional isolates from *C. plicatus* in Sa Kaeo (SK137), *H. armiger* in Khon Kaen (KK209), *H. larvatus* in Kamphaeng Phet (KP239, KP261, KP263, and KP274) were left out of the multi-locus characterization because they represented genogroups (Cp3, H1, and H3) which had already been detected. In summary, from the total of 42 isolates obtained from the bats, we narrowed the number of isolates selected for multi-locus characterization to 30 which were representative of the 7 *Bartonella* genogroups identified by *gltA* sequences (Fig 1). As noted in the main text, *Bartonella* genogroups segregated strictly among bat genera (Table A): genogroups Cp1-3 were found only in *C. plicatus*, genogroups H1-3 were found only in *Hipposideros* spp., and genogroup Tm was only found in *T. melanopogon*. However, it should be noted that the small sample size and limited species overlap among locations in our study precludes us from accurately estimating the rate at which these bat species may share *Bartonella* genogroups, especially if exchanges are rare.
**Table A. Counts of bat species captured and *Bartonella* genogroups cultured in four provinces of Thailand.** Total infected bats are recorded for each location. The counts of *Bartonella* genogroups for some species are larger than the number of infected bats because multiple cultures were isolated from some individuals, including some coinfections of multiple genogroups. Exact 95% binomial confidence intervals for prevalence were estimated using the Clopper-Pearson method.
|Province|Species|Family|Count|Positive|Prevalence (%)|95% confidence interval (%)|*gltA* genogroups (count)|
|:----------:|:----------:|:----------:|:----------:|:----------:|:----------:|:----------:|:----------:|
|Chiang Rai|*Aselliscus stoliczkanus*|Hipposideridae|3|0|`r prev.conf[2, 4]`|[`r prev.conf[2, 5]`, `r prev.conf[2, 6]`]||
||*Hipposideros fulvus*|Hipposideridae|1|1|`r prev.conf[3, 4]`|[`r prev.conf[3, 5]`, `r prev.conf[3, 6]`]|H3 (1)|
||*Rhinolophus* sp.|Rhinolphidae|1|0|`r prev.conf[1, 4]`|[`r prev.conf[1, 5]`, `r prev.conf[1, 6]`]||
||||5|1|`r loc.conf[1, 4]`|[`r loc.conf[1, 5]`, `r loc.conf[1, 6]`]||
|Kamphaeng Phet|*Hipposideros armiger*|Hipposideridae|4|0|`r prev.conf[4, 4]`|[`r prev.conf[4, 5]`, `r prev.conf[4, 6]`]||
||*Hipposideros larvatus*|Hipposideridae|23|12|`r prev.conf[5, 4]`|[`r prev.conf[5, 5]`, `r prev.conf[5, 6]`]|H1 (10), H2 (2), H3 (3)|
||*Rhinolophus* sp.|Rhinolophidae|1|0|`r prev.conf[6, 4]`|[`r prev.conf[6, 5]`, `r prev.conf[6, 6]`]||
||*Taphozous melanopogon*|Emballonuridae|1|1|`r prev.conf[7, 4]`|[`r prev.conf[7, 5]`, `r prev.conf[7, 6]`]|Tm (2)|
||||29|13|`r loc.conf[2, 4]`|[`r loc.conf[2, 5]`, `r loc.conf[2, 6]`]||
|Sa Kaeo|*Chaerephon plicatus*|Molossidae|41|17|`r prev.conf[8, 4]`|[`r prev.conf[8, 5]`, `r prev.conf[8, 6]`]|Cp1 (10), Cp2 (2), Cp3 (6)|
|Khon Kaen|*Chaerephon plicatus*|Molossidae|3|0|`r prev.conf[10, 4]`|[`r prev.conf[10, 5]`, `r prev.conf[10, 6]`]||
||*Hipposideros armiger*|Hipposideridae|14|3|`r prev.conf[9, 4]`|[`r prev.conf[9, 5]`, `r prev.conf[9, 6]`]|H1 (3), H3 (1)|
||*Hipposideros* sp.|Hipposideridae|1|0|`r prev.conf[11, 4]`|[`r prev.conf[11, 5]`, `r prev.conf[11, 6]`]||
||||18|3|`r loc.conf[4, 4]`|[`r loc.conf[4, 5]`, `r loc.conf[4, 6]`]||
|**Total**|||93|34|`r conf[, 4]`|[`r conf[, 5]`, `r conf[, 6]`]||
##### Separate gene trees for MLST loci
Maximum likelihood (ML) phylogenies generated for each locus (*ftsZ*, *gltA*, *nuoG*, *rpoB*, and ITS) separately using RAxML v8.2.10 [@Stamatakis2014] showed that most cultures clustered into the same genogroups as identified by *gltA* sequences (Figs C-G). The two exceptions were isolates KP174 and KP287a. Isolate KP174 from a *Hipposideros* sp. bat in Kamphaeng Phet was identified as genogroup H1 by all loci except *nuoG* where it clustered with genogroup H3. Isolate KP287a from a *H. larvatus* individual in Kamphaeng Phet was identified as genogroup H1 by all loci except *rpoB* where it clustered with genogroup H2.
![**Fig C. Maximum likelihood phylogeny for *ftsZ* sequences from Thai bats.** The phylogeny was created with RAxML using the GTRCAT model with 25 per site rate categories. Node support was estimated with 1000 bootstrap replicates and is indicated by the size and color of circles at each node. Tip labels for recombinant strains are colored red.](RAxML_ftsZ_bipartitions_fixed.png)
![**Fig D. Maximum likelihood phylogeny for *gltA* sequences from Thai bats.** The phylogeny was created with RAxML using the GTRCAT model with 25 per site rate categories. Node support was estimated with 1000 bootstrap replicates and is indicated by the size and color of circles at each node. Tip labels for recombinant strains are colored red.](RAxML_gltA_bipartitions_fixed.png)
![**Fig E. Maximum likelihood phylogeny for *nuoG* sequences from Thai bats.** The phylogeny was created with RAxML using the GTRCAT model with 25 per site rate categories. Node support was estimated with 1000 bootstrap replicates and is indicated by the size and color of circles at each node. Tip labels for recombinant strains are colored red.](RAxML_nuoG_bipartitions_fixed.png)
![**Fig F. Maximum likelihood phylogeny for *rpoB* sequences from Thai bats.** The phylogeny was created with RAxML using the GTRCAT model with 25 per site rate categories. Node support was estimated with 1000 bootstrap replicates and is indicated by the size and color of circles at each node. Tip labels for recombinant strains are colored red.](RAxML_rpoB_bipartitions_fixed.png)
![**Fig G. Maximum likelihood phylogeny for ITS sequences from Thai bats.** The phylogeny was created with RAxML using the GTRCAT model with 25 per site rate categories. Node support was estimated with 1000 bootstrap replicates and is indicated by the size and color of circles at each node. Tip labels for recombinant strains are colored red.](RAxML_ITS_bipartitions_fixed.png)
#### Sequence References
Sequences representing reference strains of *Bartonella* species and also previously identified *Bartonella* strains from bats or other mammal hosts were used to reconstruct the *gltA* phylogeny (Fig 1) and multi-locus phylogeny (Fig 2) in the main text. GenBank accession numbers for these sequences are listed in Table B.
**Table B. GenBank accession numbers for *ftsZ*, *gltA*, *nuoG*, *rpoB*, and ITS sequences from *Bartonella* reference strains and *Bartonella* genogroups from mammal hosts.** Blank references for some strains indicate that no sequence for that locus has been listed on GenBank. References for some *Bartonella* strains from bats and other mammals are listed.
```{r ref_table, echo=F, include=T, message=F, warning=F}
ref.table = read.csv("New_TableS2.csv")
kable(ref.table, align="c")
```
#### References