Privatization of a breeding resource by the burying beetle Nicrophorus vespilloides is associated with shifts in bacterial communities

It is still poorly understood how animal behaviour shapes bacterial communities and their evolution. We use burying beetles, Nicrophorus vespilloides, to investigate how animal behaviour impacts the assembly of bacterial communities. Burying beetles use small vertebrate carcasses as breeding resources, which they roll into a ball, smear with antimicrobial exudates and bury. Using high-throughput sequencing we characterize bacterial communities on fresh mouse carcasses, aged carcasses prepared by beetles, and aged carcasses that were manually buried. The long-standing hypothesis that burying beetles ‘clean’ the carcass from bacteria is refuted, as we found higher loads of bacterial DNA in beetle-prepared carcasses. Beetle-prepared carcasses were similar to fresh carcasses in terms of species richness and diversity. Beetle-prepared carcasses distinguish themselves from manually buried carcasses by the reduction of groups such as Proteobacteria and increase of groups such as Flavobacteriales and Clostridiales. Network analysis suggests that, despite differences in membership, network topology is similar between fresh and beetle-prepared carcasses. We then examined the bacterial communities in guts and exudates of breeding and non-breeding beetles. Breeding was associated with higher diversity and species richness. Breeding beetles exhibited several bacterial groups in common with their breeding resource, but that association is likely to disappear after breeding.

Burying beetles, Nicrophorus vespilloides, rear their offspring in small vertebrate carcasses which 59 carcasses ('beetle' treatment). Burying beetles prefer breeding in fresh carcasses over decaying ones, 112 potentially because decayed carcasses yield them lower reproductive success (Rozen et al., 2008). 113 Hence our choice of introducing beetles after 1 day of decomposition is reflective of the study 114 species' natural behaviour. The remaining 6 carcasses were manually buried (approximately 2 cm 115 deep) in the soil, to mimic the burial performed by beetles ('buried' treatment). On day 3, we removed 116 beetle pairs from the beetle treatment because carcass preparation was complete. On day 4 of the 117 experiment, we sampled the bacterial community on the beetle and buried carcasses. 118 To sample the bacterial communities on carcasses, we first removed as much soil debris as possible 119 with sterile tweezers. To maximize the amount of bacterial DNA collected from every region of the 120 carcass, we rolled the carcass in 40 ml of sterile phosphate-buffered saline (PBS) on Petri dishes, 121 using a sterile swab to release as much material as possible from the carcass into the PBS solution. 122 We pipetted the solution to a 50ml tube and pelleted the bacterial cells and debris at 3930 × g for 10 123 minutes. We discarded the supernatant and stored the pellet at -80°C until DNA extraction. 124 By analysing bacterial samples from carcasses in the fresh treatment, we were able to characterize the 125 bacterial communities present on the surface of the carcass before introduction of beetles. By 126 comparing bacterial communities on the beetle and buried carcasses, we could account for changes in 127 microbial community that were due to carcass age and burial. 128 After breeding, females were retained. Five of those females were placed with a virgin male from the 138 laboratory stock population in breeding boxes half-filled with moist compost and provided with a 139 thawed mouse carcass. Three days later, at the time of larval hatching, we collected anal exudates 140 from the breeding females, using a standard procedure we have described before (Cotter & Kilner 141 2010). Beetles readily produce a brown exudate when tapped gently on their abdomen. Exudates were 142 collected with a capillary tube and diluted in 200 µl sterile PBS. On the same day, we collected 143 exudates of four non-breeding females that had been kept in individual boxes. Breeding and non-144 breeding females were then anesthetized with CO 2 and their gut was resected. Beetle guts were placed 145 in centrifuge tubes with sterile PBS. Gut and exudate samples were stored at -80°C until DNA 146 extraction. 147

Molecular analysis 148
DNA was isolated using the FastDNA® Spin Kit for Soil (MP Bio Laboratories, Inc. Carlsbad, CA, 149 USA) and stored at -20°C until use. In order to obtain an approximate measure of bacterial abundance 150 in different carcass treatments, we performed quantitative real-time polymerase chain reaction (qRT-151 PCR) on a fragment of the 16S rRNA-encoding gene (detailed methods in Supplementary Material). 152 For library construction, we first PCR-amplified the full length bacterial 16S rRNA-encoding gene 153 (using primer pair 27F/ U1492R; Weisburg et al., 1991). We confirmed the presence of the amplicon 154 by agarose gel electrophoresis. The amplicon band was excised from the gel and the DNA was 155 extracted with the Wizard® SV Gel and PCR Clean-Up System (Promega). DNA was quantified 156 using a NanoDrop ND-1000 Spectrophotometer. To amplify the V3 region of the 16S rRNA-encoding 157 gene, a second PCR was run using 10 ng of amplicon template DNA per sample, with Illumina-158 compatible primers and PCR conditions described in Bartram et al. (2011). Amplification of the V3 159 region was confirmed by agarose gel electrophoresis and DNA was extracted from the corresponding 160 band (300 bp). High-throughput paired-end sequencing was carried out using an Illumina MiSeq 161 instrument at the DNA Sequencing Facility (Department of Biochemistry, University of Cambridge).

Procedure described in Kozich et al. (2013) and MOTHUR's Wikipedia page 165
(http://www.mothur.org/wiki/MiSeq_SOP, accessed August 2015). The quality filtering steps are 166 described in detail in the Supplementary Material. In brief, we trimmed sequences to reduce sequence 167 variation due to sequencing errors and removed sequences with many homopolymers. We aligned 168 sequences to the SILVA release 119 reference alignment and excluded those with low search scores 169 and low similarity to the template sequences. Sequences were further de-noised during pre-clustering 170 by clustering sequences with a difference of 2 or fewer nucleotides. Chimeric sequences were 171 removed using UCHIME (Edgar et al., 2011). The remaining sequences were taxonomically classified 172 by comparison against the SILVA release 119 reference database. Taxonomic assignment was made 173 at each level, given a bootstrap value greater than 80, using the Ribosomal Database Project (RDP) 174 Classifier (Wang et al., 2007). Sequences classified as Chloroplast, Mitochondria, Archaea, Eukaryota 175 or unknown at the kingdom level were removed. Uncorrected pairwise distances were calculated 176 between sequence reads. Sequences at a distance threshold of 0.03 were then clustered into 177 operational taxonomic units (OTUs), using the average neighbour algorithm (Schloss and Westcott, 178 2011). A consensus classification for each OTU was obtained. A data matrix was generated with 179 every OTU and the number of reads belonging to each sample assigned to each OTU. To control for 180 differences in the number of reads obtained per sample, we used a sub-sample of the dataset in all 181 analyses of 10760. This was chosen because the smallest number of reads in any sample was 10760, 182 obtained in a gut tissue sample of breeding beetles. 183

Statistical analysis 184
Community richness and diversity 185 Differences in observed richness and diversity between carcass treatments (measured using the 186 inverse Simpson index) were tested with ANOVA. For the analysis of community richness and 187 diversity in beetle guts and exudates, we used type of tissue (gut or exudate) and breeding condition 188 (breeding, non-breeding) as factors in an ANOVA. The inverse Simpson index was log-transformed 189 for the beetle guts and exudates data, because a Levene test indicated variance heterogeneity.

Community membership and structure 191
Non-metric multidimensional scaling (NMDS) with three dimensions was applied to Bray-Curtis 192 distance matrices, to visualize distances between samples. Differences between communities were 193 tested with PERMANOVA in R using the adonis function (vegan package, Oksanen et al., 2015). The 194 same model structure as the ANOVAs described above was used for PERMANOVA. A posteriori 195 comparisons between levels of factors were obtained by customizing the model's contrast matrix (R 196 script provided in the Supplementary Material). 197 To identify OTUs strongly associated with each treatment and combinations of treatments, we used 198 Indicator Species Analysis in R (Cáceres and Legendre, 2009). This is a standard community ecology 199 approach that takes into account both relative abundance and relative frequency of occurrence in 200 various sites (Dufrêne and Legendre, 1997). The Indicator Value is highest when all occurrences of an 201 OTU are found in a single group of sites (i.e., treatments) and when the OTU occurs in all instances of 202 that group (i.e., samples within a treatment). 203 In order to gain insight into the associations between OTUs in different carcass treatments, we used 204 network analysis. The degree of association between OTUs within each group of samples was 205 measured as Pearson's correlation coefficient. The OTU table was separated by groups ('fresh', 206 'beetle', and 'buried') and OTUs that were not present in at least 50% of the samples within each group 207 were removed. Pearson's correlation coefficients were calculated with the rcor.test function from the 208 ltm package in R (Rizopoulos, 2006). P-values were generated maintaining the false discovery rate 209 below 5% using the Benjamini-Hochberg procedure. We considered all significant Pearson 210 correlation coefficients to be network edges. Visualization of these networks was done using the R 211

Richness and α-diversity of bacterial communities
224 Across all samples we found that most OTUs occurred at very low sequence abundances. After 225 subsampling, the OTU table for soil samples contained 2653 OTUs, with just 9% of those responsible 226 for 80% of total sequence abundance. In carcass samples, after subsampling, 549 OTUs remained. 227 Just 7 of those OTUs (1.4% of the total) contributed to 80% of the total sequence abundance. 228 The highest OTU richness and diversity was observed in soil samples (Table 1). Observed richness 229 was significantly higher in buried carcasses than in fresh and beetle-prepared carcasses ( Table 1; 230 effect of treatment: F 2 =10.62, p = 0.001). Carcasses prepared by beetles and fresh carcasses showed 231 similar levels of observed richness (Table 1; Tukey post-hoc comparison: adjusted p = 0.76). 232 Bacterial taxonomic diversity, assessed as the inverse Simpson index, also differed significantly 233 between carcass treatments (effect of treatment: F 2 = 5.93, p = 0.012). Buried carcasses showed the 234 highest diversity, while beetle-prepared and fresh carcasses showed similar levels of diversity. 235 Diversity was higher in buried than fresh (adjusted p = 0.012) and beetle carcasses, although in the 236 latter case the difference was marginally non-significant (adjusted p = 0.087). Fresh and beetle-237 prepared carcasses showed similar levels of diversity (adjusted p = 0.451). Hence, despite the higher 238 concentration of bacterial DNA, indicating that beetle carcasses are a prolific ground for bacteria, beetle-prepared carcasses exhibited relatively fewer species than unprepared carcasses of the same 240 age. 241

Community composition 242
There were significant differences in bacterial community composition across all treatments. Soil 243 communities differed significantly from carcass communities (Pseudo F = 10.554, P(Perm) = 0.001; 244 Figure 2). Communities on fresh carcasses were significantly different from communities on beetle-245 prepared carcasses (Pseudo F = 9.167, Bonferroni-corrected P = 0.002). There were also significant 246 differences between beetle-prepared and control carcasses (Pseudo F = 6.071, Bonferroni-corrected P 247 = 0.002). 248 The most common bacterial phyla in soil communities were Actinobacteria and Proteobacteria 249 ( Figure 3a). Acidobacteria and Bacteroidetes were also common, yet contributed lower proportions of 250 reads. There were also 472 OTUs that could not be classified to phylum level. 251 In fresh carcasses, the majority of reads belonged to Firmicutes OTUs (Figure 3b). In one sample, 252 OTUs belonging to Proteobacteria composed approximately 50% of the observed reads, but this Alphaproteobacteria were indicator species (Table 2). subsampling, of which just 4 (1.4% of the total) contributed to more than 80% of total sequence abundance. OTU richness was overall higher in breeding beetles than in non-breeding beetles (Table  290 1; F = 13.22, p = 0.003). The lowest observed richness was found in guts of non-breeding beetles. We 291 found a significant interaction between breeding status and type of tissue (gut or exudate), largely 292 driven by a difference in richness between guts and exudates of non-breeding beetles, although 293 statistically, this difference was marginally non-significant (p = 0.080). Breeding beetles showed 294 higher diversity than non-breeding beetles (F = 12.25, p = 0.003), independent of the type of tissue 295 (guts or exudates). variability between samples, perhaps because exudates come into contact with bacteria in the beetle's 308 external cuticle. In breeding beetles, both gut and exudate samples show variability, which makes it 309 statistically difficult to differentiate between communities. As was found in beetle-prepared 310 carcasses, exudates from breeding beetles showed a high proportion of Bacteroidetes (Figure 4b). 311 Bacteroidetes were also present in three gut samples in breeding beetles. Proteobacteria and 312 Actinobacteria also comprised a small proportion of reads in breeding beetles.

325
The impact of animal behaviour on microbial communities in the environment is still poorly 326 understood. Here we show that burying beetles, which use small vertebrate carcasses as breeding 327 resources, change the bacterial communities living on those carcasses. We found that burying beetles 328 increase bacterial load when they prepare a carcass, and that there are substantial differences in 329 richness, diversity and composition of bacterial communities developing on beetle-prepared and un-330 prepared carcasses. Bacterial communities in burying beetle guts and exudates showed similarities 331 with the carcasses upon which they breed, but these resemblances were confined to the breeding 332

event, indicating transient associations between beetles and environmental bacteria. 333
Our first finding refutes the current assumption that burying beetles 'clean' and reduce the bacterial 334 load on their carcasses (e.g . Rozen et al., 2008). On the contrary, our finding reveals that carcass 335 preparation (of relatively fresh carrion) is associated with an increase in bacterial load, rather than a 336 reduction. Our results support the conclusions of a recent behavioural study examining resource use in 337 larvae of two species closely related to N. vespilloides, which suggests that the application of 338 antimicrobial exudates on carcasses has not evolved primarily as a response to microbial competition 339 for carrion (Trumbo et al., 2016). Our study did not set out to test concurrent hypothesis for the 340 evolution of antimicrobial exudates in burying beetles, yet the absence of several known pathogens in beetle-prepared carcasses lends support to the hypothesis that antimicrobial exudates serve as social 342 immunity against pathogens Cotter and Kilner, 2010). 343 Despite the increased bacterial load, we found that bacterial communities on beetle-prepared carcasses 344 showed levels of species richness and diversity as low as those associated with fresh carcasses, and 345 lower than on the carcasses we manually buried. One explanation could be that the buried carcasses 346 were colonized by soil bacteria, which are highly diverse (Figure 2a). However, beetle-prepared 347 carcasses were also buried in the soil yet exhibited less bacterial diversity. We suggest that the Gram+ bacteria. We found support for this prediction when comparing the bacterial communities 357 found on fresh carcasses with those from beetle-prepared carcasses. However, comparing beetle-358 prepared carcasses with carcasses that we buried revealed that burying beetles are capable of 359 eliminating some Gram-bacteria too. This is because the majority of the Gram-OTUs on beetle-360 prepared carcasses were Bacteroidetes, whereas on manually buried carcasses, most Gram-OTUs 361 were Proteobacteria. Proteobacteria include several insect pathogens, such as Serratia, Pseudomonas, 362 Enterobacter and Escherichia-Shigella sp. (Bulla et al., 1975), hence there is potentially strong 363 selective pressure for beetles to reduce the abundance of these bacterial groups on their breeding 364 resource. Currently, we can only speculate about what other mechanisms are in place to control 365 Proteobacteria (the pH of anal exudates is approximately 9, A.D., personal observation). 368 Alternatively, burying beetles may promote colonization of the carcass by other bacteria, which then 369 eliminate Proteobacteria. 370 Community composition on burying beetle-prepared carcasses was also changed by the addition of 371 some groups of bacteria. The presence of beetles on a carcass was associated with Clostridiales, 372 Lactobacillales and Flavobacteriales ( Table 2). The first two groups are common in the mammalian work, it would be interesting to test whether these bacteria form a mutualistic relationship with the 382 beetles by playing an active part in inter-species "bacterial warfare". Another possibility is that they 383 are commensals that benefit from the exclusion of other bacteria or from the changes in physical 384 conditions on the carcasses. 385 We used network analysis to understand more about how bacterial groups interact with each other in 386 the different carcass treatments. Our results suggest that bacterial communities on all types of 387 carcasses are poorly connected and highly modular. Given that beetle-prepared carcasses have 388 experienced a greater disturbance than control carcasses, we would expect them to show network 389 characteristics of 'stressed' communities, such as reduced correlations among species (e.g. Sun et al., 390 2013). Instead we found that the network on beetle-prepared carcasses exhibited more correlations 391 among species than the network on manually buried carcasses and was topologically more similar to 392 the network on fresh carcasses. Interactions among bacterial groups contribute to emergent properties of the microbiome in which they live (Ley et al., 2006). Our network analyses suggest that bacterial 394 communities on fresh and beetle-prepared carcasses might be more likely to exhibit emergent 395 properties than those on carcasses we buried ourselves. 396 To understand more about the role of the burying beetle in restructuring bacterial communities, we 397 examined the bacterial communities in gut and exudates of N. vespilloides females formed during and 398 after the breeding event. Communities in breeding beetles differed from (non-virgin) non-breeding 399 beetles (Figure 4), as might be expected given the changes in diet, social environment and 400 physiological state when beetles breed. One key finding was that bacterial species richness and 401 diversity was lower in the guts of non-breeding females than the guts of breeding females. This may 402 be because some aspects of the internal immune system are down-regulated during breeding (Cotter et  many of those groups, which suggests that bacterial groups common to Nicrophoridae are picked up 418 from the environment, because they feed and breed on similar species of dead vertebrate. We also found very low bacterial diversity in our non-breeding individuals, which were wild-caught but fed on 420 minced beef before being sacrificed. This supports the conclusion that at least some of the key 421 bacteria common to the gut microbiota of Nicrophorus species are environmentally acquired on their 422 breeding resource. 423

Conclusions 424
This study describes the effect of animal behaviour on bacterial communities in the environment 425 where they breed. We show that the behavioural and chemical defences of burying beetles change the 426 diversity and richness of bacterial communities on the vertebrate carcasses used as a breeding 427 resource. Our data are not consistent with the notion that these defences rid the carcass of bacteria. 428 Instead, they suggest that beetle defences induce a shift in bacterial community membership, and 429 increasing the bacterial load on the carcass in general. Our analyses further suggest that carcass 430 preparation by burying beetles involves some hitherto undescribed mechanism for elimination of 431 potential pathogenic Gram-bacteria. The challenges for future work are to investigate the 432 functionality of the bacterial community associated with breeding beetles, and whether the shift in 433 bacterial community is adaptive for burying beetles or simply a by-product of carcass preparation. 434 This work establishes the potential of the burying beetle as a system to study the manipulation and 435 assembly of microbiotas.